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Preface 


With  the  explosive  rate  of  growth  of  the  laser  tech- 
nology since  I960,  the  laser  beam  has  become  a potential 
candidate  as  a weapon  system.  In  another  development  since 
the  early  1960's  the  theories  on  optimal  estimation  and 
control  have  been  developed.  Still  another  concept — the 
conical  scan  tracking — has  been  applied  in  the  world  of 
radar  for  a number  of  decades.  The  conical  scan  and  the 
estimation  theories  are  combined  with  a digital  minicom- 
puter to  achieve  fine  pointing  and  tracking  of  a target 
with  the  laser  beam.  The  result  of  the  labors  of  this 
thesis  is  a breadboard  implementation  of  the  above  concepts 
at  the  Air  Force  Weapons  Laboratory  (AFWL) , Kirtland  Air 
Force  Base,  New  Mexico,  where  a laser  beam  "tracks"  a 
large  ball  bearing. 

The  achievement  of  the  results  presented  here  would 
have  been  impossible  without  the  assistance  and  guidance 
rendered  me  by  Major  Richard  M.  Potter,  my  thesis  advisor. 
His  devotion  to  the  duty  of  keeping  me  on  the  proper  track 
can  best  be  exemplified  by  the  patient  hours  upon  hours  of 
discussions  and  the  piles  of  chalk-dust  in  the  chalk  tray. 
Very  sincere  thanks  go  to  Lt.  Col.  James  D.  Dillow,  this 
thesis  sponsor  at  the  AFWL,  for  his  hospitality  at  the 
laboratory,  and  for  the  material  and  moral  support  that  he 


r 


■ 

I 

K 

provided  during  my  stay  there.  Additionally,  gratitude 
is  expressed  to:  Captain  Peter  A.  Maybeck  for  his  guidance 

in  the  theory  of  estimation.  Dr.  Alan  B.  Callendar  for  his 
aid  in  the  minicomputer  programming.  Captain  James  E.  Negro, 
Lt.  Stanley  R.  Robinson,  and  Dr.  Gary  B.  Lamont  for  their 
interest  and  valuable  assistance  in  this  problem. 

The  greatest  debt  I owe  is  to  my  wife  Maryann  for  her 
enormous  love  and  understanding  during  this  trying  period. 
Her  support  was  the  single  greatest  motivator  in  the 
preparation  of  this  work. 

Zdzislaw  H.  Lewantowicz 
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ABSTRACT 


An  adaptive  extended  Kalman  filter  (EKF)  is  used  in 
an  Adaptive  Laser  Optics  Techniques  (ALOT)  control  loop  to 
track  the  glint  from  a spherical  target.  A fourfold 
improvement  in  tracking  bandwidth  is  obtained  over  a con- 
ventional conical  scan  ALOT  tracking  loop  implemented  with 
the  same  hardware.  The  increase  in  tracking  bandwidth  is 
manifested  in  the  ratio  between  the  (conical  scan)  dither 
frequency  and  the  tracking  bandwidth  which  is  2.5:1  for  the 
EKF  and  10:1  for  the  conventional  analog  or  digital  glint 
tracking  schemes.  The  EKF  provides  correct  pointing  error 
estimates  over  essentially  the  entire  range  of  the  pointing 
error,  whereas  in  a conventional  conical  scan  loop  the 
demodulated  error  versus  true  pointing  error  function  is 
highly  non-linear  due  to  the  nonlinearity  of  the  glint. 

This  results  in  a demonstrated  ability  of  the  EKF  to  y0\ 
track  a larger  amplitude  disturbance  than  the  conventional  > 
conical  scan  loop.  Adaptive  schemes  are  employed  to  make 
the  estimation  algorithm  insensitive  to  fluctuations  in 
the  mean  level  of  reflected  intensity  and  to  variations  in 
the  strength  of  measurement  noise. 
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I.  Introduction 


This  report  presents  the  use  of  a Kalman  filter  in  an 
Adaptive  Laser  Optical  Techniques  (ALOT)  glint  tracking 
scheme  for  precision  pointing  of  a laser  beam  at  a target. 
The  Kalman  filter  and  control  algorithm  has  been  simulated 
on  the  CDC  6600  computer  and  implemented  in  the  ALOT  demon- 
stration laboratory  at  the  Air  Force  Weapons  Laboratory, 
Kirtland  AFB . This  estimation  and  control  algorithm 
improved  the  ALOT  control  loop  bandwidth  by  a factor  of 
four. 

The  high-energy  laser  potentially  is  an  effective 
weapon,  especially  against  airborne  targets.  Before  the 
laser  can  be  realized  as  a weapon,  however,  several  engi- 
neering problems  must  be  solved.  One  of  the  significant 
problems  is  in  the  precision  pointing  of  the  laser  beam  at 
a target.  Any  relative  high  frequency  motion  between  the 
laser  output  beam  and  the  target  will  disperse  the  energy 
over  an  area  of  the  target  larger  than  the  beam  spot  size. 
Since  target  kill  probability  is  a function  of  the  time 
average  peak  energy  density  on  the  target  surface,  it  is 
essential  to  limit  the  beam  dispersal  to  a minimum.  The 
purpose  of  the  ALOT  control  loop  is  to  eliminate  or  reduce 
this  dispersal  by  controlling  the  small  amplitude  relative 
motion  between  the  beam  and  the  target.  This  motion  can 


result  from  target  motion,  beam  jitter,  mirror  vibration, 
optical  path  bending  by  the  atmosphere,  and  any  other 
device  jitter  that  can  cause  relative  motion,  real  or 
apparent,  between  the  beam  and  the  target.  We  will  call 
this  motion  the  disturbance  process. 

Research  Goals 

The  tracking  ability  of  the  ALOT  control  loop  can  be 
discussed  in  terms  of  a ratio  between  the  dither  frequency 
and  the  bandwidth  of  the  disturbance  process  that  can  be 
controlled.  This  ratio  will  be  called  the  dither-to- 
response  ratio. 

Although  as  expressed  earlier,  the  goal  is  to  mini- 
mize the  laser  beam  energy  dispersion  over  an  area  of  the 
target,  it  is  not  possible  to  eliminate  entirely  all  of  the 
relative  motion.  Additionally,  it  is  necessary  to  induce  a 
known  amount  of  relative  motion  between  the  laser  beam  and 
the  target  before  any  information  about  the  position  of 
the  laser  beam  with  respect  to  the  target  can  be  deter- 
mined. This  induced  motion  is  called  the  dither  and 
denoted  here  as  "d" . Although  any  induced  motion  of  the 
beam  with  respect  to  the  target  is  contrary  to  the  goal  of 
minimizing  the  relative  motion,  the  dither  is  essential. 

It  is  generally  of  small  relative  amplitude  and  its  fre- 
quency is  considerably  higher  than  the  highest  disturbance 
process  frequency  to  be  controlled. 


An  accepted  measure  of  the  tracking  loop  bandwidth  is 
the  error-rejection-bandwidth.  It  is  defined  as  the  fre- 
quency at  the  -3dB  point  on  the  error  rejection  Bode  plot. 
Error  rejection  is  defined  as  20  log  [l/  (1+G)],  where  G is 
the  unity-feedback  open- loop  transfer  function.  For  this 
problem,  error  rejection  is  defined  as  20  log | x( jw) /©T ( jw) | 
where  X and  0T  are  the  pointing  error  and  the  disturbance 
process  respectively  in  the  frequency  domain. 

The  dither-to-response  ratio  will  be  used  as  a basis 
of  comparison  for  the  different  ALOT  loop  control  schemes. 

We  can  present  two  arguments  for  using  the  dither-to- 
response  ratio  as  the  basis,  assuming  that  the  dither  fre- 
quency is  the  highest  frequency  that  the  fine  pointing 
system  can  respond  to.  First,  assume  that  the  bandwidth  at 
the  disturbance  process  to  be  controlled  is  stated  as  a 
specification,  and  that  the  cost  of  the  tilt  mirror  used 
in  the  ALOT  loop  rises  with  the  required  bandwidth  of  the 
mirror.  Then  it  is  economical  to  use  a mirror  which  has  a 
bandwidth  equal  to  the  dither  frequency  and  not  greater. 
However,  the  dither  frequency,  and  in  turn  the  mirror 
bandwidth,  is  determined  by  multiplying  the  bandwidth  of 
the  disturbance  process  that  we  want  to  control  by  the 
dither-to-response  ratio.  Thus  the  bandwidth  of  the  dis- 
turbance process  and  the  dither-to-response  ratio  determine 
the  minimum  bandwidth  of  the  fine  pointing  mirror.  Sec- 
ondly, assume  that  we  have  a fine  pointing  mirror  of  a 
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given  fixed  bandwidth,  then  the  highest  disturbance  fre- 
quency that  we  can  control  is  determined  by  dividing  the 
mirror  bandwidth  by  the  dither-to-response  ratio. 

Thus,  the  dither-to-response  ratio  is  chosen  as  the 
basis  for  comparison  between  the  previously  achieved 
performances  in  the  ALOT  control  loop  and  the  performance 
realized  in  this  study. 

Analog  ALOT  Control  Scheme  (Ref  2) 

The  analog  ALOT  control  scheme  uses  the  conical  scan 
principle  employed  in  search  and  track  radars.  It  is  not  a 
new  concept  in  laser  beam  directional  control  as  evidenced 
by: 

A multidither  tracking  and  focus  control  system 
has  at  various  times  been  referred  to  by  many  names: 
'Conical  Seem'  . . . , ' Beam  Active  Tracking'  (BAT), 
or  'Adaptive  Laser  Optical  Techniques'  (ALOT) 

(Ref  6:12) . 

The  following  discussion  of  the  analog  ALOT  control  loop 
paraphrases  portions  of  Ref  2:52-57.  This  concept  uses  the 
laser  radiation  reflected  from  the  target  to  determine  the 
angular  pointing  error.  The  laser  beam  is  dithered  in  a 
small  circle  about  a nominal  pointing  direction  and  the 
pointing  error  is  deduced  from  the  resulting  changes  in  the 
intensity  of  the  reflected  radiation.  The  sensed  reflected 
radiation  from  the  target  can  be  described  as  a function  of 
the  laser  beam  pointing  direction  and  is  called  the  glint. 
The  pointing  error  in  one  axis  is  defined  as  the  angular 
deviation  x at  the  beam  from  the  center  of  the  glint, 
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boresight  error.  This  translates  to  a linear  pointing  error 
in  the  target  plane.  The  glint  reflected  from  a ball 
bearing  can  be  modeled  as  a bivariate  Gaussian  (Normal) 
function  with  a two  dimensional  pointing  error  argument. 

I(x-y)  - *max  exp(-  * V > + b (1> 

2oG 

where  I is  the  maximum  return  intensity,  a_  is  the  dis- 
max  b 

persion  of  the  Gaussian  function,  and  b is  a bias,  as  shown 
in  Figure  la.  The  x and  y arguments  are  the  two  axes  of 
the  pointing  mirror.  Equation  (1)  can  be  rewritten  as: 

2 2 

I (x,y)  = I exp  ( j)exp( + b (la) 

max  2a_^  20/ 

For  y = 0,  Equation  (1)  becomes 

2 

l(x,0)  - I exp  (-  -^)  + b (2) 

max  2 * 

In  the  single  axis  x,  the  glint  l(x,0)  is  the  resulting 
model  of  the  reflected  intensity.  This  is  generated  by 
sweeping  the  laser  beam  along  the  "equator"  of  the  ball 
bearing.  The  reflected  intensity  increases  as  the  beam 
approaches  the  "center"  of  the  ball,  and  decreases  as  the 
beam  departs  the  center.  The  single-axis  glint  then  is  a 
unimodal  function  which  is  modeled  by  a Gaussian  curve. 

The  two-axis  and  the  single-axis  glints  are  presented  in 
Figure  1;  a single-axis  experimental  glint  is  shown  in  the 
next  chapter  (Figure  7) . 
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a.  Two  Axis  Glint  Model  without  Bias  b 
(Traced  from  Ref  8:52) 


I(x) 


b.  Single  Axis  Glint  Model  with  Bias  b 
Figure  1.  Glint  Model  for  a Spherical  Target 
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The  two-axis  pointing  system  is  presented  schematically 
in  Figure  2.  Disregarding  the  y-axis  in  Figure  2,  we  can 
discuss  the  single  x-axis  control  loop.  In  the  case  of  a 
single-axis  pointing  system,  the  laser  beam  is  dithered  by 
a small  angle  about  the  nominal  pointing  angle.  The 
resulting  change  in  reflected  intensity  e(t)  is  shown  in 
Figure  3.  The  magnitude  of  e(t)  is  approximately  propor- 
tional to  the  amplitude  of  the  pointing  error  x,  provided  x 
is  less  than  the  glint  spread  aG<  This  is  shown  in  Figure 
4.  The  tested  signal  includes  e (t) , the  nominal  glint 
intensity,  a bias  due  to  background  radiation,  and  wide 
bandwidth  (white)  measurement  noise.  This  signal  is  passed 
through  the  automatic  gain  control  (AGC) . This  makes  the 
demodulation  scheme  essentially  independent  of  Imax  as 

shown  in  Figure  4.  Variations  in  I can  be  caused  by 

max  -1 

laser  power  variation,  by  a change  in  composition  of  the 
optical  transmission  medium  (e.g.,  smoke  or  haze  passing 
the  optical  path) , or  by  a change  in  the  target  range  or 
target  curvature. 

The  signal  exiting  the  AGC  is  then  demodulated  by 
multiplying  the  AGC  output  by  the  reference  dither  (oj^) 
signal.  The  output  of  the  demodulator  consists  of  the  low 
frequency  pointing  error  signal  and  harmonics  of  the  dither 
frequency  (Ref  2) . Since  only  the  pointing  error  signal  is 
desired,  the  output  of  the  demodulator  is  filtered  to  pass 
the  pointing  error  signal  and  attenuate  all  other  signals. 
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Figure  2.  Block  Diagram  of  Conical  Scan  System  (Ref  2:55) 
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Figure  3.  Effect  of  Pointing  Error  on  Reflected  Intensity  for  Dithered  Beam  (Ref  2:53) 


Figure  4.  Error  Signal  Out  of  PMT  and  Out  of  AGC  (Ref  2:67) 

The  pointing  error  signal  is  then  fed  to  the  controller, 
which  generates  mirror  commands.  This  command  signal,  in 
turn,  drives  the  mirror  to  null  the  pointing  error  in  that 
axis.  How  well  the  pointing  error  can  he  eliminated  by  the 
analog  ALOT  control  loop  depends  primarily  on  the  frequency 
content,  the  strength  of  the  disturbance  process,  and  the 
control  loop  limitations. 

The  limitations  of  the  analog  ALOT  control  loop,  as 
implemented  at  AFWL,  come  from  a number  of  sources.  Here 
the  three  most  significant  sources  are  discussed:  the 


10 


linear  model  assumption  has  a limited  region  of  validity, 
the  small  phase  margin  due  to  the  low-pass  filter  after 
the  demodulator  and  the  inherent  design  limitation  due  to 
the  lack  of  modeling  flexibility  in  this  conventional 
analog  approach. 

The  linear  model  assumption  is  reasonably  good  for 
pointing  error  |x|»  • 5oG  as  can  be  seen  in  Figure  4.  For 
pointing  errors  beyond  |x|  » *5oG  the  deviation  of  the 
curve  from  the  straight  line  increases  greatly.  When  this 
pointing  error  signal  is  passed  through  the  AGC , the 
linear  region  is  expanded  to  |x|  * 1.5aG.  This  linear 
model  works  reasonably  well  in  the  closed  loop  ALOT  for 
target-mirror  angular  disturbances  that  are  small.  How- 
ever, if  these  angular  perturbations  exceed  |x|»1.5oG 
the  large  deviation  of  the  derived  pointing  error  from  the 
true  pointing  error  causes  ALOT  loop  to  break  lock. 

The  phase  problem  due  to  the  filter  after  demodulation 
is  best  described  by  a quotation  from  Ref  2:73: 

. . . the  higher  order  the  filter,  the  sharper 
the  cutoff.  However,  higher  order  filters  cause  addi- 
tional phase  losses  and  reduce  the  system  phase  mar- 
gin. This  reduces  the  stability  margin. 

The  combined  effect  of  the  nonlinearity  and  the  phase 

problem  due  to  the  low-pass  filter  causes  a phenomenon 

called  jump  resonance.  At  the  jump  resonance  frequency  the 

ALOT  control  loop  suddenly  introduces  a catastrophic  phase 

shift  and  amplitude  loss.  This  is  shown  in  Figure  5.  The 
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Figure  5.  Measured  Closed  Loop  Frequency  Response  with  Jump 
Resonance  for  @T  * 100  sin  wt(virad)  (Ref  2:91) 

12 


I 


igyr* 


Phase  (deg) 


effect  of  the  input  amplitude  on  the  jump-resonance  fre- 
quency and  the  error  amplitude  are  shown  in  Table  I. 


Table  I 

TABULATED  JUMP  RESONANCE  FREQUENCIES  AND  ERROR  AMPLITUDES 
AS  A FUNCTION  OF  INPUT  AMPLITUDE  | ©T | 


Input  Amplitude 

|0X1 

yrad 

Jump  Resonance 
Frequency,  Hz 
(experimental) 

Error  Amplitude 
s/x  2+y2 
Virad 

30 

52 

79 

40 

48 

85 

60 

36 

45 

100 

33 

56 

150 

24 

32 

200 

21 

30 

(From  Ref  2:92) 


According  to  Ref  2:92  (Table  I),  the  response  ratio  is 
heavily  dependent  on  the  amplitude  of  the  disturbance 
process.  For  disturbance  amplitudes  A <oQ  the  closed  loop 
bandwidth  is  approximately  50  Hz  for  a dither  frequency  of 
lKHz . This  translates  to  dither-to-response  ratio  of  20:1. 
In  Ref  6:12  an  open  loop  bandwidth  of  50  Hz  was  obtained  at 
a dither  frequency  of  750  Hz  or  1 KHz,  for  a dither-to- 
response  ratio  of  15:1  and  20:1  respectively.  Although  no 
theoretical  limit  was  calculated,  control  engineers  at  the 
Air  Force  Weapons  Laboratory  (AFWL)  feel  that  the  limiting 
dither-to-response  ratio  is  approximately  10:1  for  the 
analog  system.  This  limit  is  for  a control  system  oper- 
ating in  the  approximately  linear  region  of  the  highly 
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nonlinear  error  processor.  Recent  laboratory  work  at  the 
AFWL  has  resulted  in  a dither-to-response  ratio  of  6:1. 

The  control  method  used  is  not  available  at  the  time  of 
this  writing. 

Digital  vs.  Analog  ALOT  Control 

A digital  control  scheme  was  implemented,  at  the  ALOT 
breadboard  of  the  AFWL,  which  essentially  mimicked  the 
analog  ALOT  control  loop  (Ref  7) . The  glint  tracking  idea 
was  viewed  as  a hill-climbing  task  and  a gradient  algorithm 
using  first  differences  was  implemented  using  a minicom- 
puter. The  result  was  a dither-to-response  ratio  compa- 
rable to  that  of  the  analog  ALOT  control  loop,  10:1. 
Although  the  digital  technique  did  not  significantly 
improve  the  response  ratio,  it  did  demonstrate  the  feasi- 
bility of  replacing  the  complex  array  of  analog  devices  in 
the  ALOT  control  loop  with  a digital  minicomputer.  The 
biggest  advantage  of  the  digital  approach  is  the  tremen- 
dous flexibility  afforded  by  digital  devices  and  tech- 
niques. 

In  view  of  the  limitations  in  the  analog  and  digital 
control  schemes,  another  approach  was  necessary  to  improve 
the  response  ratio.  It  was  suggested  that  a digital  sto- 
chastic estimator  may  do  a better  job  of  extracting  the 
pointing  error  information  within  the  ALOT  control  loop. 
This,  then,  is  the  beginning  point  of  this  thesis. 


Statement  of  the  Problem 


In  this  study  an  extended  Kalman  filter  which  is  an 
optimal  estimator,  in  the  minimum  mean  squared  error  (MMSE) 
sense,  is  used  to  extract  the  pointing  error  from  the  highly 
nonlinear  glint  information.  An  "optimal"  controller  is 
used  in  the  closed  loop  system  to  improve  the  dither-to- 
response  ratio  of  the  ALOT  controller.  Three  engineering 
problems  are  emphasized:  modeling  simplicity,  feasibility 

of  the  algorithm  implementation,  and  improvement  of  response 
ratio.  These  problems  are  expanded  upon  in  later  paragraphs. 

Scope . This  study  considers  the  target  motion, 
atmospheric  disturbances,  laser  beam  instabilities,  laser 
platform  vibrations,  PMT  noise,  and  computational  errors  as 
disturbance  noise  sources.  Simple  noise  models  are  used  to 
describe  the  above  disturbances,  as  other  current  studies 
strive  for  accurate  descriptions  of  these  models.  This 
study  concentrates  on  the  proper  model  of  the  nonlinear 
measurement  function,  the  treatment  of  this  nonlinear 
relationship  between  the  measured  intensity  and  the  pointing 
error,  and  reasonable  efficiency  of  the  computational  tech- 
niques. The  nonlinear  relationship  is  treated  by  reline- 
arizing about  the  current  operating  point. 

Assumptions.  In  approaching  this  problem,  some 
assumptions  are  necessary  in  the  definition  of  the  problem. 
The  detailed  noise  models  are  not  necessary  for  the 
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demonstration  of  improvement  of  the  response  ratio,  and 
simple  noise  models  are  used  for  this  purpose.  The  measured 
return  intensity  as  a function  of  the  pointing  error  is 
modeled  by  a Gaussian  (Normal)  function,  Eq  (1) . 

Goals  and  Standards.  The  goal  of  this  study  is  to 
implement  in  hardware  and  software  a digital  estimator  and 
controller  utilizing  simple  models,  but  at  the  same  time 
attain  a significant  improvement  in  the  dither- to-response 
ratio.  As  mentioned  earlier,  the  dither-to-response  ratio 
is  the  primary  standard  for  comparison  of  the  results  of 
this  research  to  previous  results  in  this  control  problem. 
Achievement  of  useful  results  in  a control  application 
depends  not  only  on  the  control  technique  applied  to  a given 
problem,  but  also  on  the  adequate  choice  of  the  mathematical 
models  used  to  describe  the  "real-world"  problem. 

The  development  of  the  problem  is  discussed  in  the 
following  sequence.  The  control  loop  model  and  the  Kalman 
filter  equations  are  presented  in  Chapter  II.  The  simu- 
lation on  the  CDC  6600  computer  and  the  implementation  on 
the  ALOT  breadboard  are  discussed  in  Chapters  III  and  IV 
respectively.  Chapter  V contains  the  results  and  Chapter 
VI  recommends  future  research  in  this  problem  area. 


II.  ALOT  System  Model  and  Kalman  Filter 


The  Kalman  filter  is  an  optimal  (in  the  MMSE  sense) 
recursive  data  processing  algorithm  which  processes  obser- 
vation data  and  arrives  at  a "best"  possible  estimate  of 
the  unknown  variables.  However,  the  "optimality"  of  this 
algorithm  is  only  as  good  as  the  "goodness"  of  the  mathe- 
matical model  describing  the  physical  ALOT  control  loop. 

For  this  reason  a great  deal  of  attention  is  given  to  the 
development  of  the  model  for  this  problem. 

ALOT  Modeling 

The  ALOT  control  loop  model  in  one  axis  is  shown  in 
Figure  6.  Perhaps  the  most  important  part  of  that  model 
is  the  glint. 

Glint.  The  glint  is  defined  as  the  intensity  of  the 
reflected  radiation  as  a function  of  the  pointing  error. 

It  depends  on  the  geometry  of  the  target  surface  and  the 
laser  beam  cross-section.  For  the  case  in  which  the  target 
is  a ball  bearing,  the  glint  is  very  well  approximated  by  a 
bivariate  Gaussian  function,  as  in  Figure  la.  (The  laser 
beam  waist  is  approximately  one  fourth  the  diameter  of  the 
ball  bearing.)  If  the  glint  is  considered  only  in  one  axis, 
i.e.,  holding  the  other  axis  argument  at  zero,  then  the 
glint  model  becomes  a single  argument  Gaussian  function  as 
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Figure  6.  ALOT  Model  and  Kalman  Filter  Control  Loop 


shown  in  Figure  lb.  The  experimental  glint  with  the  true 
Gaussian  function  superimposed,  is  shown  in  Figure  7.  This 
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Figure  7.  Glint  vs.  x-Pointing  Error 

experimental  glint  is  modeled  very  well  by  a Gaussian 
function.  The  goodness  of  this  model  can  be  best  explained 
by  the  fact  that  the  glint  is  the  result  of  a physical 
convolution  of  the  laser  beam  cross-section  and  the  reflec- 
tivity function  of  the  ball  bearing.  The  beam  intensity 
distribution  is  Gaussian.  The  spherical  shape  of  the  ball 
baaring  target  reflects  even  a uniform  density  (small  diam- 
eter) beam  in  quite  a Gaussian  manner;  this  reflection 
characteristic  is  called  the  reflectivity  function.  Then 
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when  the  ball  bearing  reflection,  that  is  closely  Gaussian, 
is  convolved  with  the  laser  beam  that  has  a Gaussian  cross- 
section,  the  resulting  glint  is  very  closely  modeled  by  a 
Gaussian  function. 

The  glint  model  also  depends  on  three  parameters  other 
than  the  pointing  errors  in  the  x and  y axis.  These  param- 
eters are  Imax  and  b,  expressed  in  volts — the  units  of  the 
photomultiplier  tube  (PMT)  sensor,  and  oG  expressed  in 
volts — the  units  of  the  input  to  the  actuators  of  the  fine 
pointing  mirror.  For  a target  of  nonvarying  dimensions 
located  at  a fixed  distance  from  the  fine  pointing  mirror 
and  illuminated  by  non-time  varying  laser  beam,  o„  and  b 

parameters  are  constant  and  measurable.  However,  I was 

max 

found  to  be  varying  in  time.  The  amplitude  of  this  vari- 
ation varied  from  day  to  day,  but  ranged  from  .2  volts  peak 
to  1.0  volts  peak.  One  sample  of  this  variation  can  be  seen 
in  Figure  8 which  represents  the  sensor  (PMT)  output  as  a 
function  of  time  for  the  case  where  the  laser  beam  is  illu- 
minating the  ball  bearing  center.  The  average  Imax  is  3.75 
v and  bias  b is  0.25  v.  The  source  of  this  I _ variation 
is  caused  by  a variation  in  the  laser  output  power.  In  a 
real  world  application  this  I variation  can  be  caused  by 
other  additional  sources,  such  as  changing  target  geometry, 
or  changing  distance  between  the  laser  tracker  and  the 
target. 


0123456789  10 


time  (sec) 

Figure  8.  Example  of  I Variation 

max 

Bias.  Continuing  clockwise  on  the  ALOT  control  loop 
in  Figure  6,  constant  bias  b is  added  to  the  glint  and  is 
shown  experimentally  in  Figure  7.  The  magnitude  of  b 
depends  on  the  strength  of  the  background  lights  and  equip- 
ment pilot  lights.  The  laser  beam  scattering  from  the 
optical  bench  optics  also  is  sensed  by  the  PMT  and  con- 
tributes to  the  magnitude  of  b. 

Measurement  Noise.  The  measurement  noise  v entering 
the  same  summer  as  bias  b in  Figure  6,  is  modeled  as  a 
white  noise  of  strength  R.  Consider  the  data  in  Figure  9 
showing  the  measurement  noise  characteristics.  Figure  9a 
shows  a very  low  level  background  noise.  The  laser  beam 
is  capped,  only  electronic  equipment  pilot  lights  are 
turned  on.  The  60  Hz  power  coming  from  these  pilot  lights 
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Figure  9.  PMT  Noise — All  Room  Lights  Off 
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Figure  9 (continued) 
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Figure  9 (continued) 
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is  evident.  The  PMT  generated  white  noise  appears  to  be 
considerably  smaller  than  .02  v and  will  be  neglected. 

Figure  9b  shows  PMT  output  with  the  laser  beam  off  target 
and  reflecting  from  a black  velvet  background.  Figure  9c 
shows  the  PMT  output  with  the  beam  positioned  at  x=0,  y=0, 
i.e.,  at  the  center  of  the  ball  bearing  target.  The  noise 
level  increased  by  a factor  of  five  over  that  in  Figure  9b. 
(Note  change  of  scale  in  the  photos.)  Figures  9d,  9e,  and 
9f  show  the  PMT  output  for  the  beam  pointing  at  x=la„, 
x=2a„,  and  x=3o^  respectively.  Although  these  photographs 
confirm  the  wide  band  of  this  measurement  noise,  it  also 
points  out  the  noise  magnitude  dependence  on  the  pointing 
error,  or  on  the  position  of  the  laser  beam  on  the  steel 
ball  target.  Additionally,  in  Figures  9d  and  9e  a low 
frequency  (approx.  0-10  Hz)  fluctuation  is  evident  and 
appears  to  be  strongly  related  to  the  low  frequency  fluctu- 
ation evident  in  Figure  7.  The  modeling  of  these  measure- 
ment noise  sources  is  deferred  to  a later  paragraph. 

Photomultiplier  Tube  (PMT) . The  next  block  in  the  ALOT 
control  loop  is  the  PMT.  This  is  a complex  non-linear 
device  that  transforms  light  energy  entering  it  into  an 
electrical  signal.  The  non-linearity  was  shown  to  be  insig- 
nificant, therefore  the  PMT  is  modeled  as  a pure  gain — Gp. 

PMT  Low  Pass  Filter.  The  next  block  in  the  ALOT  con- 
trol loop  diagram  is  the  PMT  low  pass  filter.  Its  purpose 
is  to  filter  out  high  frequencies  well  above  the  frequency 
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content  of  the  dither  and  control  loop  signals.  This  filter 
cuts  down  the  frequency  content  of  the  "white"  measurement 
noise,  and  significantly  decreasing  its  variance.  Figure 
10  shows  the  effectiveness  of  the  PMT  low  pass  filter,  and 
Figure  11  shows  the  PMT  filter  Bode  plot. 


time  (msec) 

Figure  10.  Effect  of  PMT  Low  Pass  Filter 
Upper — Raw  (A/C)  PMT  Output 
Lower — Filtered  (A/C)  PMT  Output 

A/D.  Following  the  PMT  low  pass  filter  in  the  control 
loop  is  a sampler,  or  an  analog-to-digital  converter  (A/D) . 
It  is  a 12-bit  device  with  a conversion  speed  of  approxi- 
mately 40  psec.  It  requires  an  external  clock  pulse  for 
initiation  of  the  conversion  and  has  a 0-10  volt  sample 
range. 

Kalman  Filter.  Following  the  A/D  is  the  digital  Kalman 
filter  algorithm  which  processes  the  sampled  data  to  provide 
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Figure  11.  PMT  Filter  Frequency  Response 

the  best  possible  estimates  (in  the  MMSE  sense) , ft  and 

of  the  pointing  errors  in  the  x and  y axes  respectively. 
Since  from  this  point  in  the  control  loop  both  the  x and  y 
loops  are  identical,  only  the  x channel  will  be  discussed. 
The  development  of  the  Kalman  filter  is  deferred  to  a later 
section. 

Dither  "d" . In  a conical-scan  type  modulation-demodu- 
lation scheme  a controlled  perturbation  is  applied  in  the 
form  of  a dither  to  extract  information  from  a system.  In 
this  digital  control  scheme  it  is  necessary,  as  pointed  out 
in  an  earlier  discussion,  to  induce  a controlled  motion  in 


Phase  ( ) 


each  axis  so  that  pointing  error  information  can  be 

( ) 

extracted.  This  controlled  motion  is  in  the  form  of  a 
discrete  (step)  dither  (d)  that  is  applied  simultaneously 
with  the  loop  control  input. 

Digital  Integrator,  Optimal  Controller,  and  D/A.  The 
control  signal  u(k)  = -£(k+)+d(k),  consisting  of  the  best 
estimate  of  the  pointing  error  together  with  the  dither,  is 
added  to  the  previous  mirror  pointing  direction. 

0(k+l)  = * 0m(k)  + B u (k)  (3) 

m m m m 

where  4>  =1,  and  B =1 
m m 

In  a later  section  it  will  be  shown  that  this  is  a deadbeat 
optimal  controller.  This  digital  signal  is  then  converted 
through  the  digital-to-analog  converter  (D/A) . The  D/A  is 
a 12-bit  zero  order  hold  (ZOH) , with  a range  of  -10  to  +10 
volts.  The  D/A  conversion  is  considered  instantaneous  upon 
command.  The  output  of  the  D/A  is  a piecewise  constant 
(step)  signal. 

D/A  Low  Pass  Filter.  Since  the  step  signal  ideally 
contains  the  entire  frequency  spectrum,  it  is  not  a desir- 
able signal  to  drive  the  fine  pointing  mirror.  This  mirror, 
which  will  be  discussed  in  the  following  paragraph,  has  a 
finite  bandwidth,  however,  this  finite  bandwidth  mirror  has 
resonances  at  higher  frequencies.  Thus  the  fine  pointing 
mirror  is  susceptible  to  excitation  in  its  higher  harmonics 
by  the  higher  frequencies  contained  in  the  step  function. 

To  eliminate  this  undesirable  effect  the  step  output  of  the 
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D/A  is  passed  through  the  D/A  low  pass  filter.  The  effect 
of  the  D/A  low  pass  filter  on  the  D/A  step  output  can  be 
seen  in  Figure  12.  Figure  13  shows  the  D/A  filter  Bode 
plot.  Thus  it  can  be  seen  that  the  high  frequency  exci- 
tation of  the  mirror  is  eliminated  without  attenuating  any 
dither  or  control  information.  The  bandwidth  of  this  filter 
is  350  Hz. 
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Figure  12.  Effect  of  D/A  Low  Pass  Filter 
Upper — D/A  Output  (.5v/an) 

Lower — Filtered  D/A  Output  (.2v/an) 

0 

Disturbance  Process  T.  There  are  a number  of  sources 
of  disturbances  that  may  cause  the  unwanted  movement  of  the 
laser  beam  away  from  the  desired  pointing  direction.  It 
could  be  the  vibration  of  the  target  with  respect  to  iner- 
tial space.  The  laser  optics  platform,  primarily  the  fine 
pointing  mirror,  is  susceptible  to  vibrations.  The  induced 
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Figure  13.  D/A  Filter  Frequency  Response 

relative  motion  of  the  optical  components  and  the  fine 
pointing  mirror  motion  with  respect  to  inertial  space  causes 
the  laser  beam  to  depart  from  the  desired  pointing  direction. 
Atmospheric  disturbances  can  cause  the  refraction  of  the 
beam  away  from  the  target,  and  stray  electric  signals  in  the 
mirror  actuation  system  will  cause  an  undesirable  motion  of 
the  beam.  All  of  the  above  disturbances  cause  relative 
angular  motion  of  the  mirror  pointing  direction  with  respect 
to  the  target.  Alternatively,  we  can  think  of  this  disturb- 
ance motion  as  the  motion  of  the  target  with  respect  to  the 
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laser  beam  pointing  direction.  The  important  point  is  that 
all  of  the  above  mentioned  disturbance  sources  cause 
relative  random  motion  between  the  mirror  beam  pointing 
direction  and  the  target.  This  relative  random  motion  will 
be  called  the  disturbance  process  and  symbolized  by  0T  for 
the  x mirror  tilt  axis.  The  characteristics  of  the  disturb- 
ance process  will  be  discussed  in  a later  paragraph. 

Fine  Pointing  Mirror . The  final  component  in  this 
control  loop  is  the  fine  pointing  mirror  and  the  two  actu- 
ator assembly.  The  actuators  are  piezoelectric  stacks  that 
expand  and  contract  in  response  to  the  electrical  signal 
inputs.  The  mirror  in  Figure  14  shows  the  locations  of  the 
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Figure  14.  Fine  Pointing  Mirror 
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linear  actuators,  and  Figure  15  shows  the  transformation  of 
linear  motion  to  angular  motion.  Assume  that  the  beam 


Figure  15.  Transformation  from  Electrical 
Signal  to  Linear  Displacement 

reflection  angle  is  0N  for  a neutral  x-actuator  position. 
Assusie  that  a voltage  Av  is  applied  to  the  actuator. 
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Assuming  that  the  actuator  response  is  linear  the  dis- 
placement K Ad  is  produced  where  K is  the  actuator  gain, 
in  m 

Then  the  mirror  tilt  Ay  for  small  angles  is 


Ay  = tan"1  M « M 
r r 


where  r is  the  actuator  moment  arm. 


AvK 
m 

r 


(4) 


Also,  from  Figure  12  we  have  (Snell's  law) 

20n,  » 2Gn  + AG  (5a) 

0N'  = 0n  + Ay  (5b) 


where  GN,  is  the  new  beam  reflection  angle,  and  AG  is  the 
beam  pointing  angle  change.  Then  substituting  (5b)  into 
(5a)  we  get 

AG  - 2Ay  (5c) 

However,  we  desire  the  output  Ax  in  terms  of  Av. 


Ax  = V R2+Ax2  sinAG  (5d) 

where  R is  the  distance  from  the  fine  pointing  mirror  to 
the  target.  For  small  angles  AG,  Ax  is  small  with  respect 
to  R and  Ay  can  be  approximated  by 


Ax  = RAG  * 2RAy 


^ivK 

r m 


(5e) 


The  final  relation  is 


Ax 


2RK 


m 


(5f) 


where  Ax/Av  is  the  fine  pointing  mirror  transfer  function. 

33 


l MM  I 


The  bandwidth  of  the  fine  pointing  mirror  is  high.  Bode 
plots  in  Figures  16a  and  16b  show  the  x-axis  and  y-axis 
frequency  response  of  the  tilt  mirror.  The  frequency 
responses  in  the  two  axes  are  reasonably  close  as  expected 
from  a mirror  that  is  mounted  symmetrically  on  the  actu- 
ators and  the  pivot  point.  The  significant  resonances 
appearing  above  400  Hz  may  be  the  interaction  between  the 
flexible  mirror  attachment  points  and  the  masses  of  the 
mirror  and  actuators.  Whatever  the  cause  of  these  reso- 
nances, it  is  not  desirable  to  operate  the  mirror  at  a 
frequency  above  400  Hz.  To  preclude  the  necessity  for 
modeling  of  these  resonances,  the  highest  frequency  of 
the  control  loop,  i.e.,  the  dither  rate,  was  set  at  or 
below  400  pulses  per  second.  This  corresponds  to  a basic 
dither  frequency  of  200  Hz.  Up  to  300  Hz  the  mirror  fre- 
quency response  is  essentially  flat  and  can  be  modeled  by  a 
transfer  function  consisting  of  a simple  gain — G^. 

Decoupling  of  Axes.  In  general  the  coupling  between 
the  two  axes  can  be  expected.  Two  sources  of  possible  axes 
coupling  exist:  first  is  the  mechanical  cross -coupling  of 

the  axes  due  to  the  imperfect  mirror  and  actuator  assembly, 
and  the  second  is  due  to  the  coupling  of  the  axes  in  the 
estimation  algorithm.  The  mechanical  cross-coupling  is 
weak  at  frequencies  below  500  Hz  as  observed  in  the  labo- 
ratory, and  therefore  not  modeled.  The  axes  coupling  due 
to  the  estimation  algorithm  is  significant,  especially  at 
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Figure  16.  Fine  Pointing  Mirror  Frequency  Response 
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higher  disturbance  process  (0T)  frequencies.  The  effect  of 
this  coupling  is  presented  in  a later  chapter.  However  in 
the  ALOT  model  it  is  assumed  that  at  lower  0T  frequencies 
the  cross-coupling  is  negligible  and  modeled  as  non- 
existant.  This  decoupling  is  based  on  the  premise  that  the 
dither  and  control  in  each  axis  is  independent  from  the 
other  axis.  This  is  possible  by  alternating  the  dither  and 
control  between  the  two  channels,  that  is  the  dither  and 
control  in  one  axis  are  independent  in  time  from  the  dither 
and  control  in  the  other  axis.  Because  of  this  assumed  lack 
of  axes  correlation,  the  pointing  error  estimation  and  con- 
trol can  be  accomplished  by  two  parallel  and  independent 
algorithms,  alternating  in  time. 

Kalman  Filter 

The  Kalman  filter  introduction  in  these  paragraphs  is 
paraphrased  from  Reference  3.  The  Kalman  filter  recursive 
algorithm  uses  all  available  measurements,  regardless  of 
their  precision,  to  provide  the  best  possible  estimate  of 
the  unknown  parameters.  It  does  so  in  presence  of  meas- 
urement noises  and  system  disturbances.  The  Kalman  filter 
can  be  shown  to  be  the  best  filter  of  any  conceivable  form 
if  three  restrictions  are  met.  First,  the  system  and  the 
state  measurements  must  be  linear.  Second,  the  system  and 
measurement  noises  must  be  effectively  white,  that  is,  the 
bandwidth  of  these  noises  must  be  wider  than  the  bandwidth 
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of  the  system  to  be  controlled.  Third,  these  noises  must 
be  Gaussian,  that  is  the  distribution  of  the  noise  ampli- 
tudes takes  on  the  shape  of  a normal  (Gaussian)  bell- 
shaped curve. 

Statistically  the  state  variable  can  be  described  by  a 
conditional  probability  density  function  conditioned  on  the 
measurement  history.  Since  the  Gaussian  assumption  is  made, 
this  conditional  probability  density  function  can  be  com- 
pletely described  by  the  first  and  second  moments,  or  by  its 
conditional  mean  and  variance. 

The  ALOT  control  loop  is  a nonlinear  system.  The  non- 
linearity is  in  the  relationship  between  the  pointing  error 
of  the  laser  beam  and  the  intensity  of  the  reflected  energy 
from  the  target  as  shown  in  Eq  (1) . The  dynamic  model  of 
the  ALOT  control  loop  is  linear.  The  nonlinearity  of  the 
measurement  can  be  relinearized  about  the  operating  point. 
Because  of  the  necessity  to  relinearize,  the  estimator 
takes  on  the  form  of  an  extended  Kalman  filter  (EKF) . The 
other  two  restrictions  are  fulfilled  as  the  measurement 
noise  and  the  ALOT  system  driving  disturbance  process  are 
assumed  to  be  both  Gaussian  and  white.  The  distrubance 
process  is  narrow  band,  however,  it  can  be  modeled  as  the 
output  of  a shaping  filter  driven  by  white  Gaussian  noise. 
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Simplified  ALOT  Model . Considering  only  the  x-axis. 


Figure  6 is  simplified.  The  y-axis  is  an  "image"  of  the 
x-axis.  The  PMT  and  the  D/A  low  pass  filters  are  sim- 
plified to  a gain,  as  the  bandwidths  of  these  filters  exceed 
the  bandwidth  of  interest  in  the  ALOT  control  loop.  The 
fine  pointing  mirror  and  the  PMT  are  operated  within  respec- 
tive regions  that  can  be  modeled  by  simple  gains.  The  ZOH 
will  also  be  simplified  to  a pure  gain,  at  the  discrete 
times  k,  to  keep  the  modeling  simple.  However,  it  must  be 
remembered  that  for  control  system  frequencies  approaching 
the  Nyquist  rate  limit  a significant  phase  lag  is  intro- 
duced. Combining  all  of  the  gains  into  the  fine  pointing 
mirror  gain  Gm,  the  simplified  ALOT  loop  is  presented  in 
Figure  17.  The  parameter  that  remains  to  be  discussed  is 
the  disturbance  process  0T. 

Modeling  of  the  Disturbance  Process.  The  disturbance 
process  0^  was  described  in  an  earlier  paragraph  as  the 
random  relative  motion  between  the  mirror  platform  and  tar- 
get. It  is  not  the  purpose  of  this  study  to  prepare  an 
accurate  model  of  the  rapdom  disturbance  process,  but  to 
develop  a working  model  for  the  laboratory  implementation. 

If  it  is  to  be  controlled,  this  random  disturbance 
process  must  be  band  limited,  it  is  already  assumed  to  be 
random.  This  band  limited  disturbance  process  can  be 
modeled  as  the  output  of  a shaping  filter  driven  by  white 
noise.  Only  three  questions  then  remain  to  be  answered. 


38 


Figure  17.  Simplified  Digital  ALOT  with  Kalman  Filter 


First,  what  is  the  least  order  of  shaping  filter  that  will 
be  adequate?  Second,  what  is  the  bandwidth  of  this  shaping 
filter?  And  third,  what  is  the  white  noise  strength  driving 
the  shaping  filter?  Considerable  attention  was  given  to  the 
initial  shaping  filter.  The  random  disturbance  process,  or 
the  relative  motion  between  the  tilt  mirror  and  the  target, 
consists  of  random  motions  of  the  target,  mirror,  mirror 
platform,  and  the  atmospheric  disturbances  to  the  beam 
propagation  path.  It  seemed  logical  to  think  that  a proper 
model  of  this  disturbance  process  would  be  a second  order 
shaping  filter  driven  by  white  noise  at  the  acceleration 
level  as  shown  in  Figure  18a.  Another  shaping  filter  in 
Figure  18b  is  similar  to  the  one  above,  but  has  an  addi- 
tional noise  source  at  the  velocity  level.  Finally,  in 
Figure  18c,  is  a first  order  shaping  filter  with  the  white 
noise  driving  the  velocity.  The  three  candidates  for  a 
shaping  filter  model  were  proposed  for  various  reasons. 

The  first  order  shaping  filter  is  often  an  adequate  model 
for  a particular  disturbance  process  and  is  the  simplest  to 
implement.  The  other  two  shaping  filters  are  variations  of 
disturbance  models  sometimes  used  to  describe  the  relative 
motion  between  two  airborne  vehicles.  The  shaping  filter 
that  produces  the  best  error  rejection  is  chosen  experi- 
mentally as  discussed  in  Chapter  V.  Each  of  the  above 
shaping  filters  has  been  used  in  the  Kalman  filter  algo- 
rithm, however  the  first  order  shaping  filter  was 
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. Single  Noise  Source — Second  Order  Shaping  Filter 
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Figure  18.  Possible  Shaping  Filters 

41 


experimentally  shown  to  result  in  at  least  as  good  error 
rejection  as  the  other  two.  Consequently,  the  following 
Kalman  filter  development  is  based  on  the  first  order 
shaping  filter  model.  The  bandwidth  of  the  shaping  filter 
and  the  strength  of  the  driving  white  noise  were  perturbed 
experimentally  until  best  ALOT  control  loop  response  was 
achieved. 

The  measurement  noise  v presented  earlier  from  experi- 
mental data  is  wide  band.  Therefore,  it  will  be  modeled  as 
a white  noise  of  a given  strength. 

Filter  Equations.  Combining  the  simplified  ALOT  con- 
trol loop  model  in  Figure  17  with  the  first  order  shaping 
filter  in  Figure  18c,  the  entire  two  axes  first  order  (in 
each  axis)  model  is  presented  in  Figure  19.  The  x and  y 
sub-loops  in  that  figure  are  identical  except  for  the  sample 
timing,  thus  only  one  loop  is  developed. 

The  equations  describing  a single  loop  are  given.  The 
target  motion  0,J,(t)  with  respect  to  the  mirror  base  as  seen 
in  Figure  19  is  given  by:  • 

0^(t)  = F 0^(t)  + G w'  (t) 

= -a0£(t)  + 1 w' (t)  (7) 

Where  F ■ -a,  and  G = 1 . Discretizing  this  equation: 

0^(k+l)  - *(k+l,k)  ©^(k)  + j£+1  *(k+l,T)G(T)w'(T)dx  (8) 
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:es  ALOT  Model  with  First  Order  Shaping  Filters 


FT 

Where  4>  = e , and  T is  the  sampling  period.  The  discrete 
motion  of  the  fine  pointing  mirror  with  respect  to  the 
mirror  base  is  given  by: 

Vk+1>  - Wk>  + BmVk>  <*> 

where  4>  =1,  and  B =1. 

m m 

The  discrete  input  ux(k)  contains  the  pointing  error  esti- 
mate and  the  dither  d(-l)k. 

From  Figure  19  it  is  evident  that  the  pointing  error  x is: 

x(k+l)  = 0£(k+l)  - 9m(k+l) 

and 

0m(k)  “ °T(k)  " x(k) 

Then  subtracting  Eq  (9)  from  Eq  (8)  and  substituting 

x(k+l)  = (e"aT-l)0^(k)  + x(k)-ux(k) 

+ /k+1*(k+l,T)G(x)w' (T)dT  (10)-  . 
The  relative  motion  between  the  fine  pointing  mirror  and 
the  target  can  be  arbitrarily  modeled  as  target  motion  alone 
with  respect  to  the  mirror.  The  mirror  is  being  controlled 
by  a discrete  controller.  Then  0m(k)  changes  with  respect 
to  the  mirror  base  only  at  the  sampling  times  k.  Since  it 
is  difficult  to  work  with  0,j,(k)  it  is  replaced  by  0m(k)+x(k). 
But  the  mirror  position  0m(k)  is  a random  variable  and  is  ^ 
combined  under  the  stochastic  integral  in  Eq  (10) . The 
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pointing  geometry  is  seen  in  Figure  20.  Additionally, 
the  statistics  of  the  disturbance  process  driving  target 


Figure  20.  Beam  Pointing  Angle  Geometry 
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change,  and  are  reflected  by'  a change  of  variable  from 
w' (t)  to  w(t) . Equation  (10)  then  is  rewritten: 

x (k+1)  = 4>x (k)  + Bux(k)  + /jj+1$(k+l,T)G(T)w(T)dT  (11) 


where 


-aT 


i 


and  B = -1 . 


The  measurement  equation  is 


z<k)  • W5** 


(-*  (K)  + b + v(k)  (12) 

2 / 


Where  the  bias  b is  assumed  to  be  a measurable  constant, 
and  v(k)  is  the  white  measurement  noise.  It  is  important 
to  remember  that  the  stochastic  noise  driving  the  shaping 
filter  is  no  longer  described  by  the  statistics  of  the  tar- 
get motion  with  respect  to  the  mirror  base  directly.  The 
proper  modeling  of  this  random  motion  must  include  the 
statistics  of  the  target  motion  with  respect  to  the  mirror 
base  and  subtracted  from  it  is  the  statistical  motion  of 
the  mirror  as  it  is  driven  by  known  inputs.  These  known 
inputs,  however,  are  dependent  on  the  statistics  of  the 
target  motion. 

The  filter  equations,  (Eq  13-19) , are  recursively  com- 
puted (Ref  4:187,188).  The  mean  and  the  variance  of  the 
conditional  probability  density  function  which  describe  the 
statistics  of  the  pointing  error  x(k)  are  propagated  forward 
from  time  k+  to  time  k+l~  and  then  updated  at  time  k+l+  by 
incorporating  the  measurement  at  time  k+1.  The  term  ft(k) 
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denotes  the  estimate  of  x.  The  propagation  of  the  estimate 
is  given  by 


*(k+l")  = e"aTS(k+}  + ux(k) 


(13) 


The  propagation  of  variance  P is  given  by 


P (k+l~)  = 4>2P  (k+)  + /JV G^Qdx 

- e"2aTP(k+)  + ^(l-e"2aT) 


T . 2_,2, 


(14) 


where  Q is  the  strength1 of  the  disturbance  process  0-0  . 

T m 

The  predicted  measurement  at  time  k is  given  by 

2 . A 2. 


h(*k“,5k-)  = exP 


max 


(\-  * h-  \ 
V 2a*  ) 


+ b (15) 


where  = £(k  ) and  = 0 for  single  axis.  The  linear- 
ization of  Eq  (15) , based  upon  the  predicted  pointing 
error,  is 


H<V>  * Is 


V 


Sk~ 


T- 


*k~ 


h(«k-)-b 


(16) 


The  update  (with  the  measurement)  of  the  mean  £(k+)  and  the 
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variance 


P(k+) 


is  performed  using  the  Kalman  gain  K(k) 


P(k  )H(k.-) 

K (k)  = -* (17) 

H4(k  -)P(k')  + R (k) 


where  R(k)  is  the  strength  of  the  measurement  noise  v(k). 
The  update  of  the  variance  is 


P(k+)  = 


P(k  ) 


l-K(k)H(5ik-) 


(18) 


The  update  of  the  mean  is  accomplished  at  the  very  end  of 
this  algorithm,  immediately  after  the  measurement  c (k)  is 
taken . 


*(k+)  = *(k‘)  + K (k)  c(k)  - h(ftk-) 


(19) 


The  recursive  solution  of  this  algorithm  begins  at 
the  specified  initial  conditions: 


St( 0+)  - 0 
P(0+)  = | oG2 

The  initial  guess  of  zero  pointing  error  is  better  than 
attempting  to  determine  from  the  first  measurement (s)  the 
deterministic  pointing  error.  The  algorithm  initially 
"searches"  the  target  area  for  an  intensity  measurement  that 
is  greater  than  a preset  bound,  thus  assuring  that  the  laser 
beam  is  initially  on  the  target  (on  the  glint) . The  step 
size  of  search  algorithm  is  approximately  equal  to  the 
dispersion  (loG)  of  the  glint.  The  probability  density 
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function,  that  describes  the  possible  initial  pointing  error 
based  on  the  fact  that  the  initial  measurement  5(0)  is  above 
a certain  bound,  is  nearly  a uniform  function.  For  compu- 
tational ease  it  is  assumed  that  this  uniform  density 
function  can  be  approximated  by  a Gaussian  curve.  Then  the 
pointing  error  is  unknown  but  within  certain  bounds.  If  the 
5(0)  threshold  is  set  to  coincide  with  |x|<2aG  then  this 
would  correspond  to  the  3a  point  on  the  Gaussian  probability 
density  function.  Then 

3a (0+)  = 2aG 
P(0+)  * a2 (0+)  = | oG 

is  the  initial  condition  deemed  appropriate  to  initialize 
the  algorithm. 

Pointing  Error  Sign  Estimation.  The  estimation  algo- 
rithm was  simulated  on  the  CDC  6600  computer.  During  the 
early  efforts  it  became  clear  that  the  algorithm  did  not 
properly  estimate  the  sign  of  the  pointing  error.  This  is 
the  result  of  the  single  axis  argument  x in  the  measurement 
equation,  Eq  (12),  for  y=0.  But  = ±x,  that  is,  the 
solution  has  two  roots  and  the  algorithm,  starting  from  an 
arbitrary  initial  pointing  error  x(0),  does  not  receive 
sufficient  information  to  converge  on  the  proper  sign.  In 
Eq  (19),  the  term  5 (k)-h(«k-)  is  called  the  residual.  If  the 
estimate  of  the  pointing  error  $(k~)  is  close  to  the  true 
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pointing  error,  then  the  predicted  measurement  h($k-)  is 
close  to  the  true  measurement  £(k).  As  a result,  the  resid- 
ual is  small.  Conversely,  if  the  estimate  of  the  pointing 
error  in  sign  or  in  magnitude  is  not  good,  the  residual  is 
large.  An  upper  bound  is  calculated  for  a "normal"  range 
of  the  residual  magnitude.  Then  during  the  initial  ten 
iterations  of  the  filter  the  sign  of  the  pointing  error 
estimate  St  (k  ) is  forcibly  changed  whenever  the  filtered 
residual  magnitude  exceeds  that  bound.  The  residual  magni- 
tude is  passed  through  a low  bandwidth  low  pass  filter  to 
prevent  sign  changes  for  any  single  large  residual.  This 
allows  the  sign  change  to  take  place  only  on  a consistently 
large  residual.  Once  the  filter  converges  on  the  proper 
pointing  error,  it  is  no  longer  necessary  to  keep  testing 
the  residual. 

The  extended  Kalman  filter  algorithm  described  above 
provides  us  with  the  best  possible  estimate,  in  the  MMSE 
sense,  of  the  pointing  error.  This  optimal  estimate  is,  of 
course,  confined  to  within  the  limits  of  the  goodness  of 
the  system  model  and  the  available  measurement. 

Optimal  Controller 

According  to  the  separation  theorem 

the  optimal  stochastic  controller  for  a linear 
system  driven  by  white  Gaussian  noise,  subject  to  a 

Juadratlc  cost  criterion,  consists  of  an  optimal 
inear  Kalman  filter  cascaded  with  the  optimal  feed- 
back gain  matrix  of  the  corresponding  deterministic 
optimal  control  problem  (Ref  5:16). 


This  quotation  is  summarized  by  the  key  words:  linear, 

quadratic,  Gaussian  (LQG) . 

In  the  problem  development  we  have  already  made  the 
Gaussian  assumption.  The  ALOT  control  loop  model  is  a 
linear  one,  however,  the  measurement  equation  is  highly 
non-linear.  The  extended  Kalman  filter  linearizes  this 
non-linearity  at  each  current  operating  point  and  is  hence 
treated  as  a linear  system.  A quadratic  performance  index 
"J"  is  used  to  describe  the  cost  function.  This  estimation 
and  control  problem  meets  the  requirements  of  Gaussian 
noise  and  quadratic  cost  criteria.  However,  the  relinear- 
ization about  the  current  operating  point  does  not  satisfy 
the  criterion  of  linearity.  Thus  because  of  the  lack  of  a 
better  approach  a "forced"  separation  design  is  presented, 
where  linearity  is  simply  assumed.  Thus  the  LQG  assumptions 
are  forcibly  met.  The  development  of  the  LQG  problem  is 
referred  to  Ref  5:9-18.  The  cost  function  is  written  in 
general  as: 

xT(i)X(i)x(i)  + uT(i)U(i)u(i) 

+ i xT(N+1)S  x(N+1)J  (20) 

where  X,  U,  and  S are  weighing  matrices  penalizing  the  cost 
for  deviations  in  the  state  x,  large  control  inputs  u,  and 
the  missing  of  the  terminal  state  x(N+l) , respectively.  In 
this  regulator  problem,  the  terminal  state  is  assumed  to 
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exist  at  the  infinite  time  (N=<»)  . This  implies  that  the 
terminal  condition  is  of  no  interest,  and  therefore  the 
terminal  state  weighing  matrix  S is  set  at  zero.  Without 
intermediate  steps  we  will  state  the  optimal  control  law 
for  a scalar  problem  (Ref  5:14-15): 


U* 


*( k+) ,k 


-G£(k)S(kT) 


(21) 


where 


G*(k)  = 


B K (k+1)  4> 


m 


m 


U + K (k+1) 


G*(k)  is  the  optimal  feedback  gain,  Bm  is  the  control 


4>  is  the  state  transition  con- 
m 


constant  from  Eq  (3) , 
stant  from  Eq  (3),  and  K(k+1)  is  a state  of  the  backward 
solved  Riccati  equation. 

In  this  problem  the  fine  pointing  mirror  has  a settling 
time  that  is  much  less  than  the  sample  period  used.  Thus  it 
is  reasonable  to  assume  that  infinite  control  power  is  avail- 
able as  the  mirror  can  reach  the  command  position  within  one 
sampling  period.  (This  is  a dead-beat  controller) . There- 
fore the  control  weighing  U is  equal  to  zero.  Since  $ *1, 

and  B =1,  G*(k)  becomes  unity  and  Eq  (21)  reduces  to: 

me 


u*[i(k+),k]=  ~[l]k(k+) 


(22) 


This  control  law  is  the  degenerate  case  where  U=0.  This 
result  is  intuitively  predicted  as  for  a mirror  of  a 
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negligible  mass  and  essentially  instantaneous  response, 
whenever  a best  estimate  of  the  pointing  error  £(k)  is 
available,  the  best  we  can  do  is  to  immediately  drive  the 
mirror  with  a control  in  an  opposite  direction  by  an  amount 
exactly  equal  to  the  sensed  error.  This  concept  is  clearly 
evident  in  Eq  (22)  and  is  implemented  in  the  control  law 
Eq  (3) , except  that  the  dither  is  also  included  in  the 
mirror  drive  command. 

Although  the  result  above  is  trivial,  the  fine  pointing 
mirror  will  not  always  be  as  quickly  responding,  and  the 
infinite  control  power  assumption  will  not  be  valid.  In 
that  case  the  "forced"  L-Q-G  controller  may  be  the  most 
practical  to  implement. 

The  presentation  of  the  concepts  of  application  of  the 
LQG  solution  to  the  ALOT  control  loop  has  not  assured  us  of 
improved  performance.  This  assurance  can  come  only  from  a 
simulation  on  a large  scale  data  processing  system,  such  as 
our  CDC  6600  computer  system,  or  from  an  actual  hardware 
implementation  of  the  algorithm.  Both  the  simulation  and 
the  hardware  implementation  were  accomplished.  First,  the 
simulation  techniques  are  presented  in  the  next  chapter. 
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III.  Simulation  of  Kalman  Filter  in  the  ALOT  Control  Loop 


The  simulation  was  performed  to  verify  the  theoretical 
development  of  the  ALOT  control  loop  model  and  to  gain 
insight  into  the  behavior  of  the  modeling  parameters.  The 
open  loop  simulation  was  performed  to  test  the  algorithm's 
ability  to  converge  on  the  proper  value  of  the  pointing 
error,  to  observe  the  dynamics  of  the  variance  and  Kalman 
gain,  and  to  determine  the  filter's  sensitivity  to  vari- 
ations in  dither  amplitude,  noise  strength  values,  and 
other  parameters. 

Simulation  Concepts 

A single  program  approach  is  taken  when  building  the 
simulation  program.  This  CDC  6600  computer  simulation  is 
performed  for  the  single  axis  only.  Each  significant  point 
of  simulation  techniques  is  briefly  discussed,  but  the 
results  are  deferred  until  a later  chapter. 

White  Noise  v.  The  FORTRAN  library  includes  a random 
number  generator  with  a uniform  probability  density  function 
for  numbers  between  0 and  1,  and  therefore  a mean  of  0.5. 
Twelve  random  samples  from  this  function  are  passed  through 
the  algorithm 

12 

v(k)  = a_  l [ V {12  (k-1) +i}-0 . 5J  (23) 

K i=l  u 

for  k = 1,  2,  3 . . . 
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where  v (k)  is  the  random  number,  and  oD  is  the  variance  of 
the  desired  white  noise.  This  algorithm  produces  a relative 
frequency  of  the  amplitude  that  is  a very  much  Gaussian 
distributed  white  sequence  v(k). 

Disturbance  Process  ®T  ^ . Two  approaches  are  con- 
sidered in  the  generation  of  the  disturbance  process  0T,  a 
model  of  relative  motion  between  the  fine  pointing  mirror 
and  the  target.  The  first  approach  is  to  pass  a white 
Gaussian  sequence,  similar  to  that  produced  by  Eq  (23) , 
through  a digital  first  order  filter  of  variable  bandwidth. 
The  output  of  this  filter  is  the  exponentially  time- 
correlated  disturbance  process  0T.  The  second,  and  the 
adopted  approach,  is  a generated  sinusoid  of  a given  ampli- 
tude and  increasing  frequency  in  a controllable  manner.  The 
second  approach  is  chosen  because  the  simulation  results  are 
compared  to  the  results  produced  by  a frequency  response 
analyzer  used  for  experimental  analyses.  The  frequency 
analyzer  generates  the  disturbance  process  as  a sinusoid  of 
a given  amplitude  and  varying  frequency.  It  should  be  noted 
that  the  two  approaches  would  be  equivalent,  by  super- 
position, in  a linear  system,  for  the  purpose  of  measuring 
the  error  rejection  of  the  closed  loop.  Since  in  this 
problem  a significant  non-linearity  exists  in  the  meas- 
urement of  the  pointing  error,  the  superposition  does  not 
apply.  Thus  the  two  approaches  are  only  approximately 
equivalent. 
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Frequency  Response  Algorithm.  This  algorithm  is  used 
in  generating  frequency  response  plots  presented  in  a later 
chapter.  The  algorithm  employs  the  Fourier  coefficients  in 
calculating  the  log-magnitude  of  the  response  ratio.  Given 
that  ©T(k)  is  the  controlled  sinusoidal  disturbance  signal, 
and  r(k)  the  signal  of  interest,  the  Fourier  coefficients 
are  defined: 

i N 

a (i)  = l r (k)  sin(2irf  .kT) 
k=l  l 

1 N 

bl  (i)  = Nf  I r(k)  cos  (2irf  ±kT) 
k—  1 

i N 

C1  (i)  = I ^(k)  sin(2irfikT) 

i N 

dl (i)  = 0T<k)  sin(2irfikT) 

k — 1 

where  f . is  a discrete  frequency  being  considered, 
is  the  sampling  period.  Then  the  log-magnitude  is 

a?(i)  + b?(i) 

L (i)  = 20  log  — 

m 2.2. 

c‘(i)  + d‘(i) 

where  i indicates  the  frequency  f^.  A value  N = 350  is 
chosen  for  frequencies  f^  < 80Hz. 

Operation  of  the  Simulation  Program 

The  simulation  program  is  sufficiently  flexible  to 
permit  evaluation  of  many  parameters  without  reprogramming. 
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(24a) 

(24b) 

(24c) 

(24d) 

and  T 
computed : 

(25) 


In  the  initial  phases  of  simulation , the  convergence  of  the 
Kalman  filter  in  the  ALOT  open  control  loop  is  accomplished. 
The  glint  is  stationary  for  that  simulation.  The  glint  is 
then  "disturbed"  by  driving  the  pointing  direction  with  a 
sinusoidal  0T.  In  the  final  step  the  loop  is  closed  and 
the  simulated  ALOT  control  system  performance  is  observed. 

A perturbation  analysis  is  then  performed  on  the  key  model 
parameters,  Q — the  strength  of  the  disturbance  process, 

R — the  strength  of  the  measurement  noise,  and  the  dither 
size. 

Three  different  model  versions  of  the  Kalman  filter 
are  simulated.  First,  a first-order  square-root  extended 
Kalman  filter  is  simulated,  where  the  square  root  of  the 
variance  P(k)  is  propagated  in  time.  Then  two  second  order 
filters  with  a single  noise  source  (Fig  18a)  and  with  two 
noise  sources  (Fig  18b)  are  analyzed.  Finally,  after  the 
hardware  implementation,  a first  order  filter  (Fig  18c)  is 
simulated. 

The  real  test  of  the  LQG  algorithm  is  in  its  imple- 
mentation on  the  Nova  800  minicomputer  at  the  AFWL  ALOT 
laboratory  breadboard. 


57 


IV.  Implementation  of  the  Algorithm 


The  extended-Kalman-f ilter  ALOT  control  loop  is  demon- 
strated 6n  the  AFWL  ALOT  breadboard  at  Kirtland  AFB.  The 
estimation  and  control  algorithm  is  implemented  on  a Data 
General  Corporation  Nova  800  minicomputer,  and  an  Electronic 
Associates,  Inc.  TR-48  analog  computer.  The  TR-48  is  used 
primarily  as  an  analog  interface  between  the  Nova  800  and 
both  the  fine  pointing  mirror  power  amplifiers  and  the  PMT 
sensor  output.  Using  R-C  components,  the  PMT  and  the  D/A 
filters  (Fig  6),  are  built  on  the  TR-48  patchboard.  The 
analog  computer  potentiometers  serve  as  gain  selectors  in 
the  ALOT  control  loop.  The  Nova  800  is  a 16-bit  processor 
with  an  instruction  time  (other  than  a memory-reference  or 
multiply/divide)  of  0.8  ysec.  A memory  reference 
instruction  takes  1.6  ysec,  and  a divide/multiply — 8.8  ysec. 
The  Nova  assembly  language  is  very  flexible  allowing  rela- 
tively complex  logic  operations  in  few  instructions.  The 
memory  core  is  20,000g  word  capacity,  although  only  a 
small  fraction,  less  than  7%,  of  the  core  was  occupied  by 
the  algorithm.  The  digital  control  implementation  presented 
some  special  programming  problems. 


Implementation  Problems 

In  general  the  problems  involved  in  the  implementation 
of  a complex  estimation  and  control  algorithm  are  signif- 
icant. The  complexity  and  the  size  of  the  program  gener- 
ally increase  considerably  with  the  increae  of  the  number 
of  states  in  the  control  loop  and  the  filter  model.  This 
fact  is  the  primary  motivation  for  using  simplest  possible 
models  in  the  algorithm.  Although  a second  order  filter  is 
the  highest  filter  order  implemented,  the  programming  in 
assembly  language  on  the  Nova  800  presents  a number  of 
challenges.  These  challenges  are  primarily  in  the  precision 
of  arithmetic  operations,  scaling  from  floating  point  deci- 
mal to  binary  integer  arithmetic,  round-off,  and  the  logical 
flow  of  the  digital  program.  The  methods  used  in  meeting 
those  challenges  are  presented  in  the  following  discussion. 

Precision.  The  Nova  800  has,  as  do  most  other  general 

purpose  minicomputers,  a fixed  wordlength  of  16  bits.  This 

includes  the  15  magnitude  bits  and  1 sign  bit.  The  quanti- 

15  -1 

zation  precision  of  the  least  significant  bit  is  (2  - 1)  , 

or  .0000305.  This  is  poor  by  comparison  to  the  precision  of 
large  data  processing  computers,  but  it  is  adequate  for  many 
accuracy  requirements  of  control  algorithms.  This  pre- 
cision, with  proper  scaling  (discussed  in  the  next  section) 
proved  sufficient  for  all  but  one  set  of  variables  in  this 
algorithm.  The  variance  P(k)  has  a numerical  dynamic  range 
from  10.0  to  10"6  depending  on  the  choice  of  values  of  the 
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model  parameters.  This  accuracy  is  outside  the  accuracy 
range  of  16  binary  bits,  however,  well  inside  the  accuracy 
range  of  32  bits,  i.e.,  double-word  precision.  The  arith- 
metic manipulation,  specifically  multiplication  and  division, 
involving  double  words  became  cumbersome  and  time  consuming. 

Scaling.  Due  to  a large  numerical  range  of  the  vari- 
ables, three  scale  factors,  including  a variable  scale  fac- 
tor, are  employed.  The  architecture  of  the  Nova  800  per- 
forms the  multiplication  operation  by  multiplying  the 
contents  of  two  single-words  with  the  result  becoming  a 
double  word.  Conversely,  a double-word  is  divided  by  a 
single-word  and  the  result  is  a single  word.  Then  if  a 
double-word  variance  P(k)  is  to  be  multiplied  by  another 
single-word  variable  K(k),  the  P (k)  first  has  to  be  divided 
by  an  appropriate  scale  factor,  then  multiplied  by  K(k). 

The  variable  scale  factor  is  used  in  P (k)  computation  to 
prevent  a possible  zero  intermediate  result  when  a small, 
but  significant,  double-word  number  is  divided  by  a large 
single-word  number.  A possibility  exists  of  an  intermediate 
result  being  zero,  which  in  turn  would  erroneously  render  a 
zero  final  result.  The  variable  scaling  is  accomplished  by 
shifting  left,  or  right  as  necessary,  to  perform  multipli- 
cation, or  division  respectively,  by  integer  powers  of  two. 
The  necessity  of  scaling  introduces  numerous  multiply  and 
divide  operations  in  addition  to  those  required  by  the  algo- 
rithm itself.  Specifically,  12  multiply/divide  operations 


60 


are  necessary  in  one  cycle  of  the  single  state  filter  algo- 
rithm however,  with  the  scaling/descaling  operations  the 
number  of  multiply/divide  operations  rises  to  26.  For  a 
two  state  filter  algorithm  this  number  rises  to  approx- 
imately 55.  In  the  one  state  filter  algorithm  the  single 
axis  computation  time  is  228.8  ysec  for  the  multiply/divide 
arithmetic  operations  alone.  One  cycle  computation  time 
for  the  entire  algorithm,  including  the  adaptation  schemes 
to  be  discussed  later,  was  at  approximately  750  ysec/axis 
as  compared  to  an  algorithm  sample  period  at  5 ysec/axis. 

Round-off  and  Truncation.  Truncation  error  is  intro- 
duced during  division  operations  with  finite  word  lengths. 

It  is  modeled  as  a stochastic  process  described  by  a uniform 
probability  density  function  with  a mean  of  0.5  bits.  In 
general  this  mean  is  inadvertantly  summed,  and  significant 
bias  errors  can  result.  To  eliminate  the  truncation  errors, 
the  division  result  can  be  rounded-off.  This  shifts  the 
mean  of  the  uniform  probability  density  function  to  the  zero 
point.  Thus  the  truncation  error  is  not  cumulative,  and  is 
essentially  eliminated.  Only  a few  quotients,  those  in 
sensitive  calculations,  are  rounded-off,  because  to  round- 
off all  quotients  would  result  in  a significant  computation 
time  increase. 

Exponential  Evaluation.  The  evaluation  of  the  exponen- 
tial function  poses  a peculiar  problem.  The  evaluation 
needs  to  be  accurate  as  well  as  rapid  in  computation. 
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Several  schemes  are  considered,  however,  the  scheme  chosen 
is  a table  look-up  of  exp (-2)  for  0<z£8.0.  This  provides  a 
range  for  the  pointing  error  -4lx£4,  or  for  a bivariate 
argument  -4Svx  +y  S4.  It  is  not  necessary  to  consider 
values  of  x beyond  4.0  as  the  glint  essentially  does  not 
exist  beyond  that  point,  also  the  value  of  exp(-z)=0  for 

g 

z>8.  A table  of  2 =256  entries  is  generated.  Then  each 
computed  argument  is  split  into  2 8-bit  bytes.  The  most 
significant  byte  is  then  the  relative  address  within  the 
table.  The  least  significant  byte  is  multiplied  by  the 
slope  between  the  addressed  and  adjacent  points  within  the 
look-up  table.  This  interpolation  is  then  added  to  the  base 
value  for  exp(-z)  with  accuracy  better  than  the  16  bits 
permit.  The  speed  of  this  exponential  evaluation  is  less 
than  30  ysec,  far  less  than  for  any  other  scheme  considered. 

Program  Modularity.  Due  to  the  program  complexity  it 
is  necessary  to  use  subroutines  for  essentially  all  calcu- 
lations. The  main  program  then  merely  consists  of  a number 
of  subroutine  calls.  All  subroutines  are  called  with  the 
pertinent  arguments  passed  in  a similar  manner  that  a 
FORTRAN  subroutine  call  passes  arguments.  This  is  necessary 
as  the  arguments  differ  from  one  axis  to  the  other.  The 
modularity  of  the  program  lends  itself  well  to  debugging  of 
errors  as  each  subroutine  can  be  debugged  separately. 
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Two  Axis  Program 


The  algorithm  is  coded  in  the  Nova  assembly  language. 
The  program  is  structured  to  provide  flexibility  to  the 
user. 

Program  Flexibility.  The  program  is  logically  organ- 
ized and  easy  to  use.  After  entry  of  the  few  algorithm 
parameters,  the  program  computes  all  of  the  necessary 
constants  internally  allowing  frequent  and  rapid  changes  of 
any  parameters  between  runs.  After  the  computation  of  con- 
stants, all  of  the  variables  are  initialized  as  necessary. 
The  main  program  follows  the  initialization;  it  is  composed 
mostly  of  subroutine  call  statements,  so  it  is  short  and 
easy  to  follow.  In  the  main  program,  first  the  x-axis  then 
the  y-axis  algorithms  are  computed,  then  the  main  program 
recycles  through  the  loop.  This  can  be  seen  in  the  program 
listing  in  the  Appendix. 

Monitor  Table.  Since  the  program  execution  rate  is 
much  faster  than  the  transfer  rate  of  the  available  input/ 
output  devices,  a method-  is  devised  to  record  the  time 
histories  of  variables  of  interest.  Two  tables  of  256  words 
each  are  filled  with  the  specified  variables,  thus  recording 
the  time  histories  of  the  initial  256  values  of  these  vari- 
ables. The  contents  of  these  tables  are  recorded  through  a 
hard  copy  unit  for  analyses. 

Flow  Diagram.  A functional  flow  diagram  is  presented 
in  Figure  21.  This  flow  diagram  omits  the  details  of  the 
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Figure  21.  Main  Program  Functional  Diagram 


specific  programming  language.  The  specific  details  of  the 
Nova  800  assembly  language  program  are  available  in  the 
program  listing  in  the  Appendix. 

Error  Rejection  Bode  Plots 

The  experimental  error  rejection  Bode  plot  data  are 
generated  by  an  "x-y  plotter"  attached  to  the  BAFCO  fre- 
quency analyzer.  The  BAFCO  provides  a sinusoidal  disturb- 
ance input,  ©T*  cosojt  to  the  control  loop  and  measures  the 
log-magnitude  of  the  ratio  of  the  Fourier  coefficients  of 
pointing  error  x and  disturbance  0T  at  each  discrete  fre- 
quency The  frequency  is  increased  discretely  over 

a preselected  range,  and  at  each  point  the  log-magnitude  is 
evaluated  and  plotted  (also  displayed  digitally) . 

The  effectiveness  of  the  algorithm  is  tested  in  an 
experimental  implementation.  The  results  of  the  simulation 
and  implementation  are  presented  in  the  next  chapter. 
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V.  Results 
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The  simulation  and  experiment  results  are  presented  in 
this  chapter.  The  simulation  results  establish  the  basis 
for  model  choices  and  parameter  selections  subsequently 
applied  in  the  filter  control  algorithm  on  the  ALOT  bread- 
board. 

CPC  6600  Simulation 

A significant  amount  of  simulation  was  performed  for  a 
single  axis  estimation  model.  In  the  early  stages  of  simu- 
lation, qualitative  analysis  was  performed  on  a first  order 
Kalman  filter,  to  observe  the  effects  of  different  param- 
eters on  the  quality  of  the  estimation. 

Open  Loop  Simulation.  With  the  strength  Q of  the 
disturbance  process  model  set  non-zero,  a 20  Hz  sinusoidal 
motion  0T  is  applied  to  the  simulated  fine  pointing  mirror. 
With  the  measurement  noise  strength  R set  at  a small  value, 
and  the  sampling  period  T set  at  .002  sec,  the  esti  *.tion 
is  performed  in  an  open  loop  configuration.  The  goodness 
of  the  estimation  is  observed  by  comparing  the  history  of 
the  estimated  error  to  the  true  pointing  error.  The  algo- 
rithm tracks' @T  quite  well  in  an  open  loop  configuration. 

At  this  disturbance  frequency,  the  error  in  estimation  is 
less  than  5%  of  0T  amplitude.  The  shaping  filter  is  a part 
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of  the  Kalman  filter  design,  i.e.,  different  shaping  filters 
output  different  stochastic  processes  as  models  for  0T. 

Choice  of  Shaping  Filter.  A number  of  choices  of 
shaping  filter  are  available,  as  models  for  the  disturbance 
process.  First  order  shaping  filter,  as  presented  in  Figure 
18c,  is  relatively  easy  to  implement  and  very  often  adequate 
for  a given  application.  Other  shaping  filters  such  as  the 
one  in  Figure  18a,  are  often  used  in  modeling  of  disturb- 
ances of  airborne  targets.  No  data  representing  the  real 
world  disturbance  process  is  available  for  comparison  with 
the  hypothesized  models.  The  second  order  shaping  filter 
(Figure  18a)  is  appealing  because  it  can  provide  velocity 
estimates  as  well  as  position  estimates.  It  is  intuitively 
obvious  that  a good  velocity  estimate  can  improve  a position 
estimate — a desirable  goal.  However,  the  first  order  filter 
(Figure  18c)  is  much  simpler  to  implement.  Since  the  imple- 
mentation on  the  minicomputer  deals  with  a finite  word 
length,  either  double  word  precision,  as  discussed  in  the 
previous  chapter,  or  a square-root  version  of  the  extended 
Kalman  filter  variance  has  to  be  implemented.  The  square 
root  filter  is  a Kalman  filter  version  which  propagates  the 
square  root  of  the  variance  /p (k)  (Ref  4:88-118).  This 
form  reduces  the  accuracy  requirement  at  the  expense  of  the 
algorithm  complexity.  The  square  root  version  of  the  first 
order  Kalman  filter  in  the  ALOT  closed  loop  was  simulated. 
The  filter  converges;  however,  it  estimates  poorly  for 
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other  than  low  frequencies.  The  second  order  shaping  filter 
works  reasonably  well  in  the  estimation  of  position  but 
provides  a very  poor  velocity  estimate.  The  velocity  esti- 
mate is  essentially  zero  for  all  circumstances.  A noise 
source  added  to  the  velocity  level,  thus  modifying  the 
shaping  filter  to  that  shown  in  Figure  18b,  does  not  improve 
the  velocity  estimate.  At  this  point  it  appears  that  a 
velocity  estimate  is  not  possible  as  the  velocity  state  is 
not  measurable.  Additionally,  when  the  noise  strength  Q2 
is  adjusted  to  produce  good  filter  performance,  the  value 
(including  zero)  of  the  noise  strength  does  not  alter 
the  filter  performance.  This  implies  that  only  the  inte- 
grator driven  by  white  noise  w2  is  a good  model  of  the 
stochastic  process  0T*  But  this  is  the  same  shaping  filter 
as  the  one  in  Figure  18c  with  a set  at  zero.  At  this  point 
it  can  only  be  speculated  that  the  simulation  of  the  square- 
root  version  of  the  first-order  filter  had  a programming 
error  or  other  problems. 

The  discussion  in  the  above  paragraph  presents  the 
background  information  that  led  to  the  final  implementation 
on  the  Nova  800  of  the  first-order  shaping-filter  extended 
Kalman  filter.  Double  word  precision,  instead  of  square- 
root  filtering,  is  used  to  obtain  the  necessary  accuracy  for 
the  propagation  of  the  variance  P (k) . The  above  evidence 
that  tha  two  state  shaping  filter  is  not  as  good  a repre- 
sentation of  the  real  world  disturbance  process  0T 
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characteristics  as  the  one  state  shaping  filter  dictates 
the  choice  of  the  first-order  shaping-filter  model,  such  as 
in  Figure  18c.  The  bandwidth  of  this  filter  can  be  varied 
by  choosing  the  appropriate  a and  gain. 

First  Order  Filter  Simulation.  The  first  order  shaping 
filter  in  the  ALOT  closed  loop  is  simulated,  and  the  Kalman 
filter  variance  and  gain  time-histories  are  plotted.  Also 
a log -magnitude  plot  versus  frequency  of  the  error  rejection 
is  obtained.  The  Kalman  gain  K(k)  and  filter  variance  P(k) 
are  shown  in  Figures  22  and  23  respectively.  Both  plots 
show  the  rapid  transition  from  the  initial  to  the  steady- 
state  condition.  Essentially  steady-state  values  are 
reached  by  the  filter  in  one  iteration.  This  rapid  tran- 
sition is  due  to  the  small  correlation  between  the  previous 
measurement  history  and  the  current  measurements,  as 
reflected  in  a small  value  of  the  measurement  variance  R(k) 
in  Eq  (17) . This  effect  causes  a rapid  gain  in  the  esti- 
mation confidence  as  reflected  by  this  rapid  variance 
convergence.  For  greater  values  of  R(k),  both  the  variance 
and  the  filter  gain  would  have  a slower  transient  to  steady 
state. 

In  general,  the  variance  and  gain  of  the  Kalman  filter 
for  a non-linear  system  do  not  have  a steady  state  solution, 
because  the  linearized  term  H(£  k_)  is  not  constant  as  can 
be  seen  in  Eq  (16).  Then  the  variance  P(k)  cannot  have  a 
steady  state  solution  because  H(£  k_)  appears  in  the 
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intermediate  Kalman  gain  Eq  (17) . However,  in  this  problem 

it  is  our  desire  to  zero  out  the  error  by  allowing 

u(k)  = -£(k+)  ± d (k) , where  d is  the  dither  constant.  Then 

"*qiT 

for  the  case  where  a = 0,  (e  = 1) , Eq  (13)  becomes: 

*(k+l~)  * *(k+)-£(k+)±d(k) 

= ±d(k)  = ±d. 


The  linearized  observation  H(£k~)  evaluated  at  ±d  then 
becomes  H(±dk),  where  ± implies  alternating  sign  of  the 
dither.  For  this  special  case  the  magnitude  of  H(±dk)  is 
constant,  but  the  sign  of  H(±dk)  varies  with  the  changing 
sign  of  the  dither  as  seen  in  Eq  (16) 


±dkT  1 

H(±dk)  ||h(±dk)-b 


Consequently  the  magnitude  of  the  Kalman  gain  solution  K(k) 
is  constant,  although  its  sign  changes  with  the  change  of 
the  sign  of  the  linearized  observation  parameter  H(±dk)  as 
it  is  evident  from  Eq  (17).  Since  the  magnitudes  of  H(±dk) 
and  K(k)  are  constant,  then  the  variance  propagation  P(k~) 
and  the  variance  update  P(k+)  have  steady  state  solutions, 
as  evident  in  Figures  22  and  23.  Filter  stability  and 
dynamics  are  important  in  the  system  analysis,  but  the 
effectiveness  of  the  Kalman  filter  in  the  ALOT  closed  loop 
needs  to  be  measured. 
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Error  Rejection.  The  measure  of  "goodness"  of  the  ALOT 


control  loop  is  the  ratio  of  the  dither  frequency  to  the 
error-rejection  bandwidth.  The  error  rejection  has  been 
obtained  using  the  method  presented  in  Eqs.  (24)  and  (25) . 
The  simulated  control  loop  error-rejection  curve  is  plotted 
in  Figure  24.  The  initial  slope  is  20  dB  per  decade  indi- 
cating a first  order  response.  The  -3dB  point  occurs  at 
41  Hz — indicating  that  the  error-rejection  bandwidth  is 
41  Hz.  The  dither  frequency  of  100  Hz  (200  samples  per 
second)  is  used  in  this  simulation.  The  dither  to  response 
ratio  is  2.44:1.  For  frequencies  above  40  Hz,  the  Bode 
plot  indicates  higher  order  dynamics  that  are  probably 
induced  by  a combination  of  filter  dynamics  and  non-line- 
arity effects.  No  analysis  has  been  performed  to  determine 
the  source  of  this  phenomenon.  Phase  angle  is  not  computed 
for  this  simulation.  Additionally,  this  one  axis  simulation 
of  the  ALOT  loop  is  not  tuned  for  best  performance.  The 
parameters  I , oQ,  Q,  R,  and  dither  size  are  approximately 
the  same  as  the  parameters  in  the  experimental  control  loop. 
In  the  simulation  these  parameters  are  perturbed  in  the 
small  to  maximize  the  ALOT  control  loop  bandwidth.  These 
parameters  are  not  perturbed  in  the  large  to  globally 
maximize  these  parameters  for  a maximum  bandwidth  perform- 
ance as  this  is  not  within  the  scope  of  this  feasibility 
study. 
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Figure  24.  Simulated  Error  Rejection 

ALOT  Laboratory  Results 

The  algorithm  first  is  implemented  in  single  axis  using 
a pipe  for  a target.  The  pipe  is  held  in  a vise  with  the 
pipe  axis  perpendicular  to  the  laser  beam  plane  of  move- 
ment. The  glint  from  the  laser  beam  sweeping  across  the 
pipe  is  modeled  very  well  by  a Gaussian  function. 
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Single  Axis  Results.  The  results  of  the  single  axis 


Kalman  filter  implementation  are  discussed  in  this  section. 
The  error  rejection  plot  in  Figure  25  shows  that  the  band- 
width (-3dB  point)  is  62  Hz.  This  bandwidth  translates  into 
a dither-to-response  ratio  of  1.61:1.  Considerable  time  was 
spent  "tuning"  the  filter  and  control  loop  to  obtain  this 
response.  It  must  be  remembered  that  the  theoretical  limit 
is  at  best  a dither-to-response  ratio  of  1:1  due  to  Shanon's 
sampling  theorem.  The  dither-to-response  ratio  of  1.61:1  is 
better  than  the  2.44:1  obtained  in  simulation.  However,  it 
must  be  remembered  that  simulation  algorithm  was  not  nearly 
as  well  tuned  as  the  experimental  filter.  The  important 
point  to  be  observed  is  that  the  shape  of  the  experimental 
log-magnitude  plot  (Figure  25)  compares  quite  favorably  with 
the  shape  of  the  simulated  plot  (Figure  24) . 

The  time  history  of  the  variance  is  shown  in  Figure  26. 
This  plot  shows  that  the  variance  reaches  steady  state  with- 
in two  iterations,  very  much  like  the  simulated  variance  in 
Figure  23.  The  experimental  gain  time  history  in  Figure  27 
shows  essentially  steady  state  Kalman  gain  magnitude.  This 
is  an  important  point.  In  this  algorithm  the  discrete-time 
Riccati  equation  is  solved  for  each  iteration  to  provide  the 
Kalman  gain.  However,  if  the  Kalman  gain  is  essentially 
constant,  then  it  can  be  approximated  by  a constant,  and 
thus  reduce  the  algorithm  computation  time  by  approximately 
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50%  as  well  as  eliminate  the  wordlength  problem  in  the 
variance  computation. 

The  "goodness"  of  the  single  axis  estimation  and  con- 
trol algorithm  is  presented  in  the  following  series  of 
photographs  of  signals  recorded  on  a storage  oscilloscope. 
Figure  28a  shows  the  amplitude  and  frequency  of  the  disturb- 


ance sinusoid  0T  = .2  sin(2ir  lOt)  volts  (smooth  curve). 
Superimposed  on  top  of  it  is  the  feedback  signal  0^  with  the 
discrete  dither  added.  The  dither  frequency  is  100  Hz 
(sampling  period  T = .005  sec).  The  mean  of  0m  appears  not 
to  deviate  from  0T.  This  indicates  that  the  laser  beam  is 
tracking  the  center  of  the  steel  pipe  (top  of  the  glint) 
despite  the  10  Hz  disturbance.  This  fact  is  indeed  evident 

I 

I in  Figure  28b.  The  feedback  signal  0m  is  shown  in  the 

bottom  of  the  photograph;  the  change  in  both  the  vertical 

and  horizontal  scales  for  Figure  28a  to  Figure  28b  should  be 

noted.  In  the  top  of  the  photograph  is  the  output  of  the 

PMT  filter,  or  the  measurement  sampled  by  the  A/D  converter. 

I is  set  equal  to  8 volts.  In  this  photo  the  PMT  filter 
max 

output  does  not  indicate  that  the  laser  beam  is  precisely 
at  the  top  of  the  glint  which  would  be  represented  by  8 
volts,  but  it  indicates  the  effect  of  the  dither.  The 
dither  causes  the  laser  beam  to  "jump"  on  either  side  of  the 
glint  peak.  The  important  point  to  notice,  however,  is  that 
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Figure  28.  Tracking  Photographs  (Non-Adaptive) 
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uniform,  confirming  the  observations  in  Figure  28a  that 
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This  means  that  essentially  the  only  pointing  error  present 
is  due  to  the  dither.  Of  course  this  does  not  hold  true  for 
higher  frequencies  of  ©T.  Figure  28c  shows  that  the  feed- 
back signal  ©m  does  not  follow  ©T  at  feT  = 30Hz.  This  is 
also  evident  in  Figure  28d  where  the  control  applied  at  the 
beginning  of  each  sampling  period  does  not  evenly  "jump" 
about  the  mean.  The  three  steep  ramps  between  sampling 
times  are  created  by  the  motion  ©T  causing  the  beam  to  move 
away  from  the  top  of  the  glint  and  reduce  the  magnitude  of 
reflected  intensity.  This  change  in  pointing  error,  of 
course,  cannot  be  sensed  until  the  next  sampling  time,  at 
which  time  this  deviation  is  sensed  and  the  ©m  is  changed  to 
restore  the  beam  to  the  top  of  the  glint.  At  fQT  = 60Hz  the 
estimation  errors  are  significant,  and  the  magnitude  of 
movement  of  the  laser  beam  away  from  the  top  of  the  glint 
between  the  sampling  times  is  quite  large.  This  is  approxi- 
mately the  frequency  at  which  the  error  rejection  curve 
reaches  the  -3dB  log-magnitude.  At  this  frequency,  despite 
the  above  estimation  and  control  errors,  the  control  loop 
is  tracking  and  is  stable. 

It  should  be  pointed  out  that  the  initial  sign  esti- 
mation, in  this  single  axis  implementation,  is  performed 
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through  the  technique  discussed  in  Chapter  II.  That  is, 
during  the  initial  ten  iterations  the  pointing  error  esti- 
mate sign  is  forced  to  the  opposite  sign  if  the  residual 
history  exceeds  a certain  bound.  This  works  quite  well  as 
no  hesitations  in  convergence  on  the  glint  were  observed. 
Additionally  the  algorithm  is  restarted  when  it  ceases 
tracking.  To  help  determine  the  ALOT  control  lock-on  the 
glint,  the  residuals  and  the  measurements  are  continually 
monitored.  Both  are  filtered  with  low  pass  digital  filters 
of  low  bandwidth  to  extract  the  time  averaged  values.  The 
residual  RES (k)  is  first  rectified  then  filtered  to  provide 
SRES(k)  (filtered  |RES(k)|) 

SRES(k)  = — — T | RES  (k)  | 

1-.986z"A 

The  bandwidth  of  the  low  pass  filter  is  approximately  1Hz, 
so  that  only  slow  variations  are  not  attenuated.  Similarly 
the  measurements  £ (k)  are  filtered 

SZET(k)  = — y C (k) 

1-.977z“a 

where  SZET(k)  is  the  filtered  measurement.  The  bandwidth  of 
this  filter  is  approximately  1.5Hz.  The  break- lock  monitor 
restarts  the  estimation  and  control  algorithm  whenever 
SRES (k)  exceeds  a predetermined  upper  bound  indicating 
excessively  large  estimation  errors,  or  whenever  SZET(k) 
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exceeds  a predetermined  lower  bound,  thus  indicating  that 
the  laser  beam  is  not  on  a measurable  part  of  the  glint. 

The  low-pass  filter  constants  and  the  bounds  were  chosen 
with  a trial  and  error  method. 

After  the  discussion  of  the  single  axis  implementation, 
it  is  necessary  to  examine  the  results  of  the  two  axis 
problem. 

Two-Axes  Implementation . The  modeling  assumptions  for 
a spherical  target  made  in  Chapter  II  allow  a simplified 
approach  to  this  problem.  This  simplified  approach  is 
evident  in  Eqs.  (13) -(19)  where  one  set  of  estimation 
equations  for  the  x-axis  are  presented.  An  identical  set 
of  estimation  equations  are  used  for  the  y-axis  except  that 
in  Eqs.  (13),  (17),  (18),  and  (19)  the  variable  x is 
replaced  by  the  variable  y,  also  Eq  (16)  becomes  a partial 
derivative  of  with  respect  to  y and  evaluated  at 

y*^-.  This  approach  permits  the  decoupling  of  the  disturb- 
ance processes  0T  in  the  x and  y axes,  which  in  turn  permits 
to  model  the  control  problem  with  two  parallel  control 
loops.  It  must  be  pointed  out  that  the  pointing  error  esti- 
mates are  coupled  through  the  measurement  Eq  (15) . However, 
this  coupling  was  not  modeled  in  the  covariance  propagation 
equations.  This  coupling  is  shown  in  the  following  dis- 
cussion to  be  weak  for  low  frequency  disturbance  0T,  and  it 
becomes  significant  with  increasing  0T  frequency.  It  is 
suspected,  but  not  shown,  that  this  coupling  is  responsible 


for  the  reduced  dither-to-response  ratio  for  the  two  axes 
as  compared  to  the  single  axis  where  the  coupling  problem 
does  not  exist.  Each  of  these  loops  with  its  own  estimator 
and  controller  would  time  share  the  computational  capacity 
of  the  Nova  800  minicomputer.  The  dither  scheme  imple- 
mented moved  the  beam  in  a square  box  pattern  in  a counter- 
clockwise direction.  One  complete  path  about  this  square  is 
® defined  as  one  cycle  of  the  dither.  In  the  two-axis  imple- 
mentation the  dither  frequency  is  chosen  as  100  Hz — the  same 
dither  frequency  as  for  the  single  axis  case.  Since  no 
simulation  of  the  two  axis  control  was  performed,  the  com- 
parisons will  be  made  with  the  single  axis  performance  and 
with  the  performance  of  other  tracking  schemes. 

The  first  two-axis  data  presented  are  the  error 
rejection  Bode  plots  in  Figure  29  for  the  x and  y axes. 

These  error  rejection  plots  were  generated  by  applying  0T 
to  one  of  the  axes  without  disturbing  the  other  one.  Figure 
29  shows  the  error-rejection  bandwidth  is  34  Hz  and  36  Hz 
for  the  x and  y axes  respectively.  The  average  dither-to- 
response  ratio  for  the  two  axes  is  then  2.86:1.  No  clear 
explanation  was  yet  found  for  the  almost  50%  reduction  in 
bandwidth  from  the  single  axis  results.  Perhaps  the  cross 
axis  uncoupling  assumption  in  the  Kalman  filter  equations 
is  breaking  down  at  disturbance  frequencies  around  30  Hz 
and  above. 
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Figure  30  (continued) 


similar  to  that  for  the  single  axis  estimation  was 
attempted,  but  its  success  rate  was  extremely  poor.  As  a 
last  resort  all  sign  forcing  schemes  were  abandoned;  the 
result  was  a total  success.  It  is  possible  that  the  algo- 
rithm may  not  succeed  in  the  estimation  of  the  sign  each 
time,  however,  its  success  rate  is  excellent.  It  is  not 
possible  to  perceive  algorithm  restarts  visually  due  to 
failure  in  correct  sign  estimation,  either  because  there 
were  no  failures  or  because  the  few  failures  were  masked  by 
the  rapidity  of  the  algorithm  restart.  It  is  possible  to 
set  up  a monitor  within  the  algorithm  to  observe  the  numbers 
of  sign  estimation  failures,  but  this  was  not  done. 

ALOT  Model  Parameter  Perturbations . In  the  ALOT  closed 
loop  set  up,  the  significant  model  parameters  are  perturbed 
and  their  effects  are  plotted.  To  uniformly  measure  the 
effect  of  perturbations  of  the  parameters  the  disturbance 
(0T)  frequency  is  held  constant  at  20  Hz  and  the  log-magni- 
tude from  the  BAFCO  is  noted.  An  increase  in  the  log-magni- 
tude indicates  a decrease  in  performance  of  the  ALOT  control 
loop,  and  conversely,  Figures  31  through  33  show  the  effects 
of  changing  the  values  of  the  model  parameters  from  the 
nominal.  Specifically  the  plots  show  perturbations  of  dis- 
turbance process  strength  Q,  measurement  noise  strength  R, 
and  dither  step  size.  It  is  evident  from  these  plots  that 
individual  perturbations  of  these  parameters  produce  local 
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minima  that  correspond  to  highest  error  rejection  bandwidth. 
No  global  perturbation  analysis  was  accomplished  on  the 
model  parameters. 

Target  Search  and  Acquisition.  A search  and  acqui- 
sition scheme  was  implemented  to  permit  rapid  target  search 
and  lock  on.  The  search  window  was  set  at  the  limits  of 
the  field  of  view  imposed  by  the  ALOT  optics.  The  window 
size  is  1.2  mrad  by  1.2  mrad,  where  for  comparison  the 
target  glint  dispersion  (aQ)  is  42  yrad.  A spiral  search 
pattern  was  chosen  for  its  implementation  simplicity.  The 
initial  search  is  begun  at  the  window  center  and  spirally 
diverged  outward  until  the  outer  boundary  is  reached  at 
which  time  the  search  is  reinitiated  at  the  center.  If  the 
return  intensity  at  any  of  the  discrete  search  points  is 
above  a preset  threshold,  the  algorithm  control  is  handed- 
off  to  the  estimation  and  control  algorithm.  If  for  any 
reason  the  control  algorithm  breaks  lock  on  the  target, 
control  is  returned  to  the  search  algorithm.  At  this  point 
the  spiral  search  is  initiated  at  the  last  mirror  pointing 
position.  Again  if  a boundary  is  reached  before  target  is 
acquired,  the  search  is  reinitiated  at  the  origin.  This 
scheme  is  not  optimal,  but  it  was  found  to  be  effective  for 
this  demonstration. 

Algorithm  Sensitivity.  The  performance  of  the  algo- 
rithm was  shown  to  be  dependent  on  the  correct  choices  of  Q, 
R,  and  dither  step  size.  Incorrect  selection  of  other 
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parameters,  such  as  I or  an,  also  degrades  the  perform- 
ance  of  the  algorithm.  This  is  especially  true  if  one  or 
more  of  these  other  parameters  are  time  varying.  At  least 
one  parameter  was  observed  to  vary  with  time  significantly 
in  the  ALOT  laboratory.  The  peak  return  intensity,  I , 
contained  a random  variation  of  amplitude  approximately 
equal  to  15%  of  I itself.  This  variation  is  discussed 
in  Chapter  II.  The  source  of  this  variation  appears  to  be 
the  laser  itself.  Vibrations  of  the  laser  cavity  optics 
induced  by  the  building  vibrations  appear  to  cause  the 
laser  power  output  to  fluctuate  a corresponding  amount. 

The  I fluctuations  are  strongest  whenever  the  laboratory 

IQ&X 

air  conditioning  system  is  operating.  Whenever  the  I x 
fluctuations  are  at  a minimum,  the  best  error  rejection  is 
achieved  in  single  or  in  dual  axis.  The  estimation  and 
control  algorithm  is  most  sensitive  to  Imax  fluctuations. 

The  algorithm  is  also  sensitive  to  the  assumed  strength 
of  the  wideband  measurement  noise  v.  Figure  34  shows  the 
effect  of  the  strong  wide  band  noise  added  to  the  glint 
when  the  algorithm  model  did  not  account  for  it.  Despite 
small  disturbance  frequency,  fgT  = 10Hz,  the  algorithm  is 
seen  in  Figure  34  to  commit  significant  estimation  errors 
as  compared  to  errors  in  evident  in  Figure  30b  where  the 
measurement  noise  v is  much  smaller.  The  variance  R of  the 
measurement  noise  v is  set  at  .01  volts  while  for  this 
large  amount  of  noise  R = .20  volts^  is  more  appropriate. 
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Figure  34.  Effect  of  Strong  Measurement  Noise  v 

The  strong  adverse  influence  of  Imax  variations  and 

changes  of  measurement  noise  v can  be  attenuated  in  a number 

of  ways.  In  this  experiment  l__v  and  R (strength  of  v) 

max 

were  adapted  for  using  conventional  methods  as  presented  in 
the  next  paragraph. 

Adaptive  Extended  Kalman  Filter 

Statistical  methods  are  available  for  parameter  esti- 
mation, however,  other  methods  were  applied  in  this  adap- 
tation scheme. 

Measurement  Noise  Adaptation . When  the  strength  Q of 
the  disturbance  process  is  zero,  and  it  is  bo  modeled,  then 
it  is  easy  to  show  that  the  statistics  of  the  residual  (RES) 
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are  precisely  the  statistics  of  the  measurement  noise  v when 
the  Kalman  filter  reaches  steady  state  operation. 


E{RES(k) }=0 
E{RES  (k)RES(l)  } 


0 for  l^k 
R for  l=k 


It  can  also  be  shown  that  when  the  disturbance  process 
strength  is  small,  and  the  algorithm  provides  "good"  esti- 
mates for  Qj*0,  the  following  approximation  is  valid: 

E{RES(k) }«=0 

f - 0 for  l^k 

E{RES (k) RES (1) } < 

L*  R for  l=k 

Then  under  the  assumption  of  ergodicity  the  variance  R can 

be  computed,  in  an  approximation,  as  an  average  over  a 

o 

finite  time  interval  of  RES  (k) . The  averaging  is  accom- 

2 

plished  by  passing  RES  (k)  through  a low  pass  filter  of 

approximately  1Hz  bandwidth  thus  allowing  convergence  on 

2 

only  the  slowly  varying  RES  (k)  signal.  The  output  of  the 
low  pass  filter  is  the  estimate  of  the  measurement  noise 

A 

variance  R' (k) . By  trial  and  error  it  became  evident  that 

A 

it  is  necessary  to  add  a bias  RQ  to  R' (k)  to  achieve  good 
error  rejection.  This  scheme  is  presented  in  Figure  35. 
The  gain  m of  this  low  pass  filter  is  a proportionality 
constant.  The  effect  of  adaptation  for  R can  be  seen  in 
Figure  36  where  the  estimation  errors,  and  in  turn  control 
errors,  are  significantly  smaller  than  in  Figure  34  where 


o 


the  estimation  errors  are  larger.  Despite  its  crudeness, 
the  adaptation  scheme  proved  to  be  effective  as  in  the  above 
comparison  of  the  Figures  36  and  34  and  as  will  be  shown 
later  in  improvement  of  the  dither-to-response  ratio.  The 
time  history  of  R(k)  for  one  particular  experiment  can  be 
seen  in  Figure  37  where  the  adaptation  algorithm  converges 

A o 

on  R(k)»  .17  volts  . 

Adaptation  for  jnax.  When  the  estimation  of  the 
pointing  error  is  "good" , then  given  the  pointing  error 
estimate  k(k  ) and  the  corresponding  measurement  C (k) , I , 
can  be  reconstructed.  It  is  known  that 


«<*>  - exp|-(x2+y2)/2oG2 


+ v (k) 


(26) 


Where  x and  y are  the  true  time  varying  pointing  errors. 

Solving  for  I : 

max 


Wk> 


C (k)  - v(k) 

2 — 2 2T 

exp[ - (x^+y z ) /2o(r 


(27) 


But: 


2 2 2]  k (#.-,$.-) 

.xp|-(x2+y2,/23G2J  - 


Then: 


U.-tk-UcU)  I’  (k-l)v(k) 

(28) 


Where  (k)  is  the  approximation  of  Im_„(k). 
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To  eliminate  the  high  frequency  variations  of  Imax  00  due 
to  the  white  sequence  v(k)  in  the  second  part  of  the  above 
equation,  I1Jiax(k)  is  passed  through  a low  pass  filter, 
where  S is  the  filter  gain. 


Wk>  = 


. 023S 


1-.977Z 


-1  ^ax^) 


(29) 


The  bandwidth  of  this  low  pass  filter  is  approximately 
1.5Hz.  Thus,  Im!,„(k)  responds  only  to  low  frequency 

luaX 

A 

variations  of  I'  (k) . Then  Im„„(k)  is  used  in  the  filter 

algorithm.  With  this  algorithm  implemented,  the  laser  beam 

output  was  varied  from  .5  volts  to  9.0  volts  at  the  PMT 

filter  output.  This  Imax  variation  produced  essentially 

no  variation  in  the  error  rejection  at  f0T  = 10Hz.  The 

adaptation  scheme  for  Imax  made  the  algorithm  essentially 

insensitive  to  I variations.  This  insensitivity  to  I 

max  max 

increased  the  error  rejection  bandwidth  of  the  ALOT  control 

loop  from  35  Hz  non-adaptive  to  40  Hz  with  R and  I _ 

max 

adaptations.  This  can  be  seen  in  Figure  38.  The  improved 
ALOT  error-rejection  bandwidth  translates  to  a dither- to- 
response  ratio  of  2.5:1  for  the  adaptive,  compared  to  2.9:1 
for  the  non-adaptive  estimation  scheme.  This  dither-to- 
response  ratio  is  the  best  documented  for  the  two  axis  ALOT 
estimation  and  control  in  this  research  effort.  The  time 

A A 

history  of  Imax00  i*  seen  in  Figure  39  where  Iffiax  (k) 
appears  to  be  following  slowly  varying  Imltv  (k) . The  Kalman 
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gain  and  the  variance  time  histories  for  the  adaptive 
extended  Kalman  filter  are  presented  in  Figures  40  and  41 
respectively.  The  abnormality  in  each  of  these  figures 

during  the  first  sampling  period  is  due  to  the  filter 

* 

starting  delay  of  one  iteration  to  allow  taking  of  one 
measurement  necessary  to  begin  Imax  adaptation. 

The  discussion  and  the  data  presented  in  this  chapter 
in  no  way  indicate  that  the  best  dither-to-response  ratio 
has  been  achieved.  However,  it  is  shown  that  significant 
improvement  in  the  dither-to-response  ratio  is  feasible 
through  the  use  of  a digital  adaptive  extended  Kalman 
filter. 
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VI . Conclusions  and  Recommendations 


Conclusions 

This  study  shows  that  an  optimal  estimation  scheme  can 
improve  significantly  the  response  bandwidth  of  the  ALOT 
closed  loop  tracker.  The  estimation  scheme  is  an  extended 
Kalman  filter  algorithm  implemented  on  a digital  mini- 
computer . 

The  bivariate  Gaussian  (Normal)  function  is  an  excel- 
lent model  for  the  glint  from  a spherical  target  and  a 
Gaussian  beam.  However,  the  glint  need  not  be  Gaussian  for 
the  algorithm  to  work  well.  Any  unimodal  and  continuously 
differentiable  function  that  is  a good  model  for  the  glint 
will  suffice.  This  algorithm  will  track  any  glint  that  can 
be  modeled  by  a unimodal  and  continuously  differentiable 
function.  A simple  (first  order  pure  integrator)  shaping 
filter  driven  by  white  Gaussian  noise  is  an  adequate  model 
of  the  disturbance  process  that  the  ALOT  tracking  loop 
attenuates.  This  disturbance  attenuation  permits  the 
pointing  of  the  laser  beam  at  the  top  of  the  glint.  Even 
though  the  glint  is  highly  nonlinear  for  large  pointing 
errors,  the  algorithm  can  estimate  these  pointing  errors 
close  to  the  glint  edges.  Thus  the  estimation  scheme  is  less 
dependent  on  the  amplitude  of  the  disturbance  process  than  a 
typical  conical  scan  demodulation  scheme. 
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The  non-adaptive  Kalman  filter  can  be  implemented  using 
constant  gains,  and  constant  predicted  measurements.  This 
is  due  to  the  fact  that  the  linearization  is  accomplished 
about  a fixed  magnitude  of  the  pointing  error.  This  would 
eliminate  the  need  for  the  time  consuming  re-linearization 
about  the  nominal  pointing  error  and  real  time  solution  of 
the  discrete-time  Riccati  equation. 

The  ALOT  estimation  and  control  algorithm  is  sensitive 
to  variations  of  some  parameters  that  were  assumed  to  be 
constant.  This  variation  reduces  the  bandwidth  of  the  ALOT 
control  loop.  An  unsophisticated  adaptation  scheme  used  to 
estimate  I and  R makes  the  control  loop  essentially 

IucLX 

insensitive  to  these  parameters. 

The  pointing  error  has  four  possible  sign  solutions  in 
the  two  arguments  for  any  one  measurement.  The  algorithm 
requires  no  additional  information  in  correctly  identifying 
the  signs  of  the  two  pointing  error  arguments.  This  behav- 
ior is  not  totally  explained  and  requires  further  research. 

This  research  also  shows  conclusively  that  a general 
purpose  minicomputer  can  be  utilized  for  implementation  of 
this  estimation  and  control  algorithm.  The  problems  of 
computational  errors  due  to  finite  word-length,  truncation, 
and  scaling  can  be  overcome  adequately.  The  most  important 
advantage  of  digital  control  is  the  tremendous  flexibility, 

compared  to  analog  implementation,  it  afforded  in  the 

< 

initial  programming  and  in  the  subsequent  program  changes. 
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Recommendations 


I 

t 

I 
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Further  research  is  recommended  in  the  use  of  an  adap- 
tive extended  Kalman  filter  in  an  ALOT  glint  tracking 
scheme.  Host  significant  should  be  an  analysis  of  full 
estimation  scheme  for  a given  intensity  measurement 
I(x,y,ar,I  ,b).  In  a general  "real  world"  situation 
these  parameters  are  either  variables  or  constants  not 
known  a priori.  The  goal  would  be  then  to  estimate  param- 
eters essential  to  meet  the  control  specifications.  These 
parameters  should  be  modeled  as  state  variables  in  a sto- 
chastic estimation  scheme. 

Another  area  of  suggested  research  is  in  the  correct 
modeling  of  the  cross-correlation  of  the  pointing  error 
between  the  two  axes.  The  ALOT  control  loop  performance 
can  possibly  be  improved  by  correct  modeling  of  this  cross- 
correlation. 

Offset  tracking  is  a concept  where  a portion  of  the 
unimodal  glint  other  than  the  peak  is  used  as  a tracking 
point.  This  concept  may  have  "real  world"  applications 
and  could  be  explored  in  future  research. 

The  two  axes  estimation  should  be  simulated  to  explore 
the  effect  of  cross-axis  coupling.  A proper  stochastic 
model  that  includes  the  cross  coupling  parameters  may 
improve  the  two  axis  performance  of  this  estimation  scheme. 


o 
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The  adaptation  scheme  for  the  measurement  noise  vari- 
ance R is  crude,  but  effective  to  a limited  degree.  Other 
adaptation  schemes  should  be  explored  to  improve  the  R 
estimation. 

A study  of  the  possible  targets  should  be  performed  to 
determine  the  adequacy  of  this  simplified  modeling  form. 
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Appendix 


This  appendix  presents  the  assembly  language  program 
for  the  Data  General-Nova  800  minicomputer.  The  program 
is  the  implementation  of  the  two-axes  estimation  and 
control  algorithm  of  the  digital  conical  scan  tracker. 


r 


! 
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•TITL  AXIS 2 ;THIS  PROGRAM  ESTIMATES 

;THE  LASER  BEAM  POINTING  ERROR 
7 IN  TWO  AXIS.  THIS  ERROR  IS  USED 
7 AS  A CONTROL  INPUT  IN  A FEED- 
; BACK  LOOP  TO  ZERO  OUT  THAT  ERROR. 


.EXTN 

DEBUG 

.ENT 

QQ 

.EXTU 

. ZREL 

; ALL  CONSTANTS  AND  VARIABLES 

ALPH: 

0 

;»2*PI*F 

IMAXO: 

10240. 

; INITIAL  IMAX 

SIGNT: 

1400. 

; GLINT  SIGMA  **«****##**#« 

T: 

4. 

; SAMPLING  PERIOD  *«**#*#*# 

RO: 

100. 

; INITIAL  R *««#******##*** 

CTO: 

1. 

.•INITIAL  COUNT  #***##«*### 

D: 

3200. 

; DITL  FACTOR  «####**##*#*# 

Q: 

200. 

; NOISE  STRENGTH  (DIST)  *** 

GI: 

0 

; IMAX  FLTR  GAIN  (1-GIX) 

GII: 

2048. 

7 FLTR  GAIN  MULTPLR 

GIX: 

1800. 

; IMAX  FLTR  TIME  CONST 

GR: 

0 

; R FLTR  GAIN  (1-GRX) 

GRI: 

2048. 

;FLTR  GAIN  MLTPLR 

GRX: 

1800. 

;R  FLTR  TIME  CONSTANT 

GS: 

0 

; SZET  FILTER  GAIN 

GSX: 

1024. 

;SZET  FILTER  TIME  CONSTANT 

TABC : 

0 

; STORE  TABLE  COUNTER 

TABCO: 

-301 

.•INITIAL  TABLE  COUNTER 

SF11: 

2048. 

; SCALE  FACTORS 

SF3 : 

8. 

• N 

9 

SF7 : 

16. 

N 

9 

SF15: 

32768. 

• M 
9 

SFOUT : 

10. 

;D/A  SCALE  FACTOR 

BLB: 

600. 

; BREAK  LOCK  BOUND 

PO: 

32. 

;P(0)  MSW 

P01: 

0 

; LSW 

TRES: 

1600. 

; RES  FLTR  TIME  CONSTANT 

GRES: 

0 

,-GAIN  OF  RES  FLTR 

RPROP : 

18432. 

•REST  PROPORTIONALITY  CONST 

G: 

0 

; GQG  MSW 

Gl: 

0 

7 LSW 

A: 

0 

; EXP (-ALPH*T » 

B: 

0 

?A**2 

DITL: 

400. 

; DITHER  SIZE 

SHFT: 

11. 

.•SHIFT  COUNT* iR 

SIGS: 

0 

;SIGNT**2 

SIGS2: 

0 

;2*SIGS 

UX: 

0 

; X-AXIS  CONTROL 

DAX: 

0 

7 X-AXIS  OUTPUT  REGISTER 

UY: 

0 

; Y-AXIS  CONTROL 

DAY: 

0 

7 Y-AXIS  OUTPUT  REGISTER 
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M. 


IMAX: 

0 

R: 

0 

BOUND: 

8000. 

MSK: 

177400 

MSK1: 

177740 

MSK2 : 

177000 

MSK4 : 

140000 

MSKB : 

174000 

XM: 

0 

XP: 

0 

YM: 

0 

YP: 

0 

PXM: 

0 

PXM1: 

0 

PXP: 

0 

PXP1 : 

0 

PYM: 

0 

PYM1 : 

0 

PYP: 

0 

PYP1 : 

0 

KX: 

0 

KY: 

0 

AK: 

0 

HX: 

0 

HY: 

0 

HHX: 

0 

HHY: 

0 

AHH: 

0 

ZITA: 

0 

SZET: 

0 

XRES: 

0 

SXRES : 

0 

AXRES : 

0 

YRES: 

0 

SYRES: 

0 

AYRES: 

0 

HP: 

0 

HP1: 

0 

TEM: 

0 

TEMl : 

0 

TEM2: 

0 

TEM3 : 

0 

TEM10: 

18000. 

TEM11: 

0 

TEMl 2 : 

0 

TEMl 3: 

0 

TEM14 : 

0 

COUNX: 

0 

COUNY: 

0 

ZERO: 

0 

SGN: 

0 

; CURRENT  I MAX 
; CURRENT  R 
; . 9*IMAX 
; BYTE  MASK 
; SHIFT  MASK 
; SHIFT  MASK 

; EXPONENT  ARGUMENT  MASK 
.•BOUNDARY  MASK 
; X-ERROR  AT  T- 
; X-ERROR  AT  T+ 

; Y-ERROR  AT  T- 
; Y -ERROR  AT  T+ 

; X- VARIANCE  AT  T-  MSW 
; LSW 

; X-VARIANCE  AT  T+  MSW 
; LSW 

; Y-VARIANCE  AT  T-  MSW 
; LSW 

; Y-VARIANCE  AT  T+  MSW 
; LSW 

; X-FILTER  GAIN 
; Y-FILTER  GAIN 
; ABS (KX  OR  KY) 

;H (X.Y)  IN  X 
; H (X.  Y)  IN  Y 
.‘LINEARIZED  H ABOUT  X 
; LINEARIZED  H ABOUT  Y 
; ABS (HHX  OR  HHY) 

; MEASUREMENT 
; SUM  OF  ZETAS 
; RESIDUAL  IN  X 
; SUM  OF  X RESIDUALS 
; ABS (XRES) 

; RESIDUAL  IN  Y 
; SUM  OF  Y RESIDUALS 
; ABS (YRES) 

;TEMP  STORE  FOR  H*P  MSW 
; LSW 

; SUBROUTINE  RTN  ADDR  STORE 
; TEMPORARY  STORAGE 

. » 

i 

} 

; 

> " 


.•CONVERGENCE  TIME  COUNTER 
; " 

.•DUMMY  LOCATION 
;SIGN  FLAG 
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ARG: 

0 

; EXPONENTIAL  ARGUMENT  (X**2+Y**2) 

EXPT: 

0 

; VALUE 

OF  EXPONENT 

TAB: 

0 

; TEMPORARY  STORAGE 

ZER: 

0 

; VALUE 

OF  ZERO 

ONE: 

1. 

; VALUE 

OF  ONE 

BIAS: 

820. 

; I«IMAX*EXP (- (X**2+Y**2) /SIGS2) ) +BI 

PTSIX: 

1050. 

; SRES  FACTOR 

DITIV: 

0 

; INITIAL  DITHER  SIZE 

SLPE: 

1092. 

; BOUND 

SLOPE 

MIN: 

750. 

; BOUND 

BIAS 

NO: 

1. 

; INITIAL  # OF  STEPS 

N: 

0 

; CURRENT  # OF  STEPS 

NC: 

0 

; STEP  COUNTER 

TEST: 

2000. 

; SEARCH 

THRESHOLD 

STP: 

2000. 

;STEP  SIZE 

MPP: 

MP 

IW: 

IV 

MP35M: 

MP35 

• 

9 

.NREL 

; COMPUTE  CONSTANTS  **************** 

. BLK 

20 

QQ: 

LDA 

1, ALPH 

;LOAD  ALPH 

LDA 

2,T 

; LOAD  T 

SUBO 

0,0 

;CLR  AC0 

MUL 

• 

9 

STA 

1 , TEM1 

; STORE  ARGMNT 

JSR 

@EXP 

;A=EXP(-ALPH*T) 

TEM1 

;EXP  ARGUMENT 

QQ1: 

A 

; STORE 

RSLT  IN  A 

.BLK 

10 

LOA 

1,A 

; LOAD  A 

MOV 

1,2 

; DUPLCT  A 

SUBO 

0,0 

;CLR  AC0 

MUL 

• 

9 

LDA 

2 , SF11 

; LOAD  SF11 

DIV 

• 

9 

MOVZL 

0,0 

; ROUND  OFF 

SUBZL* 

0,2,SZC 

• n 

INC 

1,1 

• tt 

9 

QQ2: 

STA 

1,B 

; STORE  B 

LDA 

2, ALPH 

; LOAD  ALPH 

MOVL# 

2,2,SZB 

;IS  ALPH-07 

JMP 

QQ3 

;NO,  JMP  QQ3 

LDA 

1,T 

; YES, LOAD  T 

LDA 

2,SF11 

; LOAD  SF11 

SUBO 

0,0 

;CLR  AC0 

MUL 

; 

LDA 

2,Q 

; LOAD  Q 

MUL 

> 

JMP 

QQ33 

; JMP  AND  STORE  G 

QQ3: 

LDA 

1, SF11 

; LOAD  SF11 

o 

QQ33: 

QQ4: 

QQ5: 

0 

QQ6: 

QQ7: 


LDA 

2,B 

SUB 

2,1 

LDA 

2 , SF11 

SUBO 

0,0 

MUL 

LDA 

2, ALPH 

DIV 

MOVZR 

1,1 

LDA 

2,Q 

SUBO 

0,0 

MUL 

STA 

0,G 

STA 

1,G1 

JMP 

QQ4 

. BLK 

20 

LDA 

1, SIGNT 

MOV 

1,2 

SUBO 

0,0 

MUL 

LDA 

2,SF11 

DIV 

STA 

1, SIGS 

MOVZL 

1,1 

STA 

1,SIGS2 

LDA 

1 ,GIX 

LDA 

2 , SF11 

MOV 

2,3 

SUB 

1,2 

LDA 

1 ,GII 

SUBO 

0,0 

MUL 

LDA 

2 , SF11 

DIV 

STA 

1 ,GI 

LDA 

1 ,GRX 

SUB 

1,3 

MOV 

3,1 

LDA 

2,  GRI 

SUBO 

0,0 

MUL 

LDA 

2,SF11 

DIV 

STA 

1 ,GR 

LDA 

1,GSX 

LDA 

2,SF11 

SUB 

1,2 

STA 

2 ,GS 

LDA 

1,SF11 

LDA 

2, TRES 

SUB 

2,1 

STA 

1,GRES 

LOAD  B 
1-B 

LOAD  SF11 
CLR  ACO 

LOAD  ALPH 

DIV  BY  2 
LOAD  Q 
CLR  ACO 

STORE  G MSW 
LSW 


; LOAD  SIGNT 
; DUPLCT  SIGNT 
;CLR  ACO 

•LOAD  SF11 

; STORE  SIGS 
;MULT  BY  2 
; STORE  SIGS* *2 
; LOAD  GIX 
;LOAD  SF11 
; DUPLCT  SF11  IN  AC 3 
; 1-GIX 
; LOAD  GII 
;CLR  ACO 

•LOAD  SF11 

; STORE  GI 
; LOAD  GRX 
; 1-GRX 

;MOV  AC 3 TO  AC1 
;LOAD  GRI 
;CLR  ACO 

; LOAD  SF11 

; 

; STORE  GR 
; LOAD  GSX 
; LOAD  SF1I 
> 1-GSX-GS 
; STORE  GS 
; LOAD  SF11 
; LOAD  TRES 
> 1-TRES 
; STORE  GRES 


O 
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LDA 

1,SIGNT  ;LOAD  SIGNT 

LDA 

2 , D ; LOAD  D 

SUBO 

0,0  ; CLR  ACO 

MUL 

# 

LDA 

2, IMAXO  ; LOAD  IMAXO 

DIV 

t 

STA 

1 , DITIV  ; STORE  DITIV 

SUBO 

0,0  ;CLR  ACO 

STA 

0,DAX  ; RESET  X-AXIS 

STA 

0,DAY  ; RESET  Y-AXIS 

IV: 

SUBO 

0,0  ; INITIALIZE  THE  VARIABLES 

STA 

o,ux 

STA 

0 , UY 

STA 

0,XM 

STA 

0,XP 

STA 

0,  YM 

STA 

0 , YP  ; 

IV1: 

STA 

0 , ZETA  ; 

IV2 : 

LDA 

1, IMAXO  ; 

STA 

1 , IMAX  ; 

STA 

1,SZET  ; 

LDA 

1,R0  ; 

STA 

1#  R 

LDA 

0,P0 

LDA 

lfPOl 

STA 

0,PXP 

STA 

0,PYP 

IV3: 

STA 

1,PXP1  ; 

STA 

1,PYP1  ; 

LDA 

1 ,CT0 

STA 

1 , COUNX  ; 

STA 

1 , COUNY  ; 

LDA 

1, DITIV  ; 

STA 

1 , DITL  ; 

STA 

lfUX 

STA 

1 ,UY 

LDA 

1 ,TABC0  ; 

STA 

1 , TABC  ; 

IV4: 

LDA 

1,MP7 

STA 

1 ,MP8  j 

LDA 

1, IMAXO  ; 

MOVZR 

1,1  ; 

STA 

1 , SXRES  ; 

STA 

1 , SYRES  ; 

JMP 

MP9+4 

IV5  j 

JSR 

0SRCH  ;FIND  THE  GLINT 

IV6: 

JMP 

MP 

. BLK 
. ZREL 

20 

P20 : 

MP20 
• NREL 
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MP: 


MP4: 


JSR 

PXM 

PXM1 

PYM 

PYMl 

JSR 

XP 

UX 

XM 

PXP 

PXP1 

PXM 

PXM1 


@TABLS 


ePROP 


MAIN  PROGRAM  ********************* 
STORE  264  SPECFD  VLS 


PROPAGATE  X-AXIS 
SUBROUTINE  ARGUMENTS 


.BLK 

20 

MP1 : 

JSR 

§ARGM 

COMPUTE  ARGUMENT 

JSR 

0EXP 

COMPUTE  EXPONENTIAL 

ARG 

ARGUMENTS 

EXPT 

.BLK 

10 

MP2 : 

JSR 

@HHHH 

COMPUTE  H,HH 

XM 

ARGUMENTS 

HHX 

HX 

.BLK 

10 

MP3: 

JSR 

0UPDT 

COMPUTE  K , P 

PXM 

- 

ARGUMENTS 

PXM1 

KX 

PXP 

PXP1 

HHX 

. BLK 

JSR 

XRES 

AXRES 

KX 

XM 

XP 

SXRES 

HX 


10 

0CLCK 


SAMPLE  AND  COMPUTE: 
XP,  XRES,  SXRES 


.BLK 

10 

JSR 

@SZETA 

; COMPUTE  SUM  OF  ZETAS 

MP5 : 

JMP 

MP8 

; BYPASS  TWO  DUMMIES 

MP6 : 

JMP 

MP10-2 

• 

9 

MP7 : 

DSZ 

COUNX 

• 

9 

MP8: 

DSZ 

COUNX 

; DECR . CONVERGENCE  COUNTER 

JMP 

MP9 

; COUNT  NOT-0,  CONTINUE 

LDA 

1,MP6 

; COUNT-0,  SET  UP  BYPASS 

STA 

1 ,MP8 

. » 
, 

JMP 

MP10 

; CONTINUE 

O 
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MP9 : 

JSR 

9BOUN 

JSR 

SXRES 

XP 

0TEN 

LDA 

1 ,MP27 

STA 

1 , MP  2 8 

JMP 

IV5 

• BLK 

5 

LDA 

2, DITL 

STA 

2 , UX 

JMP 

MPll 

MP10 : 

LDA 

1,XP 

LDA 

2,  DITL 

SUB 

1,2 

STA 

2 ,UX 

MPll: 

LDA 

1,  DAX 

ADD 

2,1 

STA 

1 , DAX 

MP12 : 

JSR 

ZER 

DAX 

@OUTPT 

.BLK 

10 

MP13 : 

LDA 

1 , COUNX 

MOV 

1,1, SZR 

JMP 

MP20 

MP14 : 

JSR 

HX 

@IEST 

JSR 
AX  RES 
RPROP 

0REST 

.BLK 

10 

MP20 : 

JSR 

YP 

UY 

YM 

PYP 

PYP1 

PYM 

PYM1 

9PROP 

.BLK 

20 

MP21: 

JSR 

9ARGM 

JSR 

ARG 

EXPT 

@EXP 

.BLK 

10 

MP22 : 

JSR 

YM 

HHY 

HY 

9HHh.  % 

.BLK 

10 

MP23: 

JSR 

9UPDT 

COMPUTE  BOUND 
IF  COUNT  NOT=0 , FORCE 
ON  XP,  OR  YP 

RESTART  Y-AXIS  BYPASS 


LOAD  DITL 
STORE  UX=DITL 
JMP  MPll 
LOAD  XP 
LOAD  DITL 
UX=-XP+DITL 
STORE  UX 
LOAD  DAX 
DAX=DAX+UX 
STORE  DAX 
OUTPUT  X CHANNEL 


LOAD  COUNT 
IS  COUNT=0? 
NO,  JMP  MP20 
ESTIMATE  IMAX 

ESTIMATE  R 


PROPAGATE  Y-AXIS 
ARGUMENTS 


COMPUTE  ARGUMENT 
COMPUTE  EXPONENTIAL 


; COMPUTE  H,HH 


5 COMPUTE  K,P 
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SIGN 


i 


i 


1 


PYM 

PYMl 

KY 

PYP 

PYP1 

HHY 

. BLK 

10 

MP24 : 

JSR 

0CLCK 

YRES 

AYRES 

KY 

YM 

YP 

SYRES 
HY 
• BLK 

10 

JSR 

6SZETA 

MP25 : 

JMP 

MP28 

MP26 : 

JMP 

MP30-2 

MP27 : 

DSZ 

COUNY 

MP28 : 

DSZ 

COUNY 

JMP 

MP29 

LDA 

1 ,MP26 

STA 

1 ,MP28 

JMP 

MP30 

MP29: 

JSR 

0BOUN 

JSR 

@TEN 

SYRES 

YP 

.BLK 

10 

LDA 

2, DITL 

STA 

2 ,UY 

JMP 

MP31 

MP30 : 

LDA 

1 , YP 

LDA 

2, DITL 

SUB  J 

1,2 

STA 

2 ,UY 

MP31 : 

LDA 

1,DAY 

ADD 

2,1 

STA 

1 ,DAY 

MP32: 

JSR 

0OUTPT 

ONE 

DAY 

.BLK 

10 

MP33 : 

LDA 

1, COUNY 

MOV 

1,1, SZR 

JMP 

MP35 

MP34 : 

JSR 

0IEST 

HY 

JSR 

@REST 

AYRES 

SAMPLE  AND  COMPUTE 
YP,  YRES , AND  SYRES 


COMPUTE  SZET 
; BYPASS  TWO  DUMMIES 

• H 
t 

• II 

;DECR  CONVER  COUNTER 
; COUNT  NOT=0 , CONTINUE 
;COUNT=0 , SET  UP  BYPASS 

• II 

; YES , JMP  MP30 
; COMPUTE  SZET 
;NO,  FORCE  SIGN 


LOAD  DITL 
STORE  UY=DITL 

LOAD  YP 
LOAD  DITL 
UY—YP+DITL 
STORE  UY 
LOAD  DAY 
DAY=DAY+UY 
STORE  DAY 
OUTPUT  Y-AXIS 


LOAD  COUNT 
IS  COUNT- 0? 
NO,  JMP  MP35 
ESTIMATE  IMAX 

ESTIMATE  R 


RPROP 

' 

. BLK 

10 

MP35: 

JSR 

0DITS 

COMPUTE  DITL  AND  NEG  IT 

LDA 

1 , COUNX 

LOAD  X-COUKTER 

MOV 

1,1, SNR 

IS  X-COUNTER=0? 

JSR 

0 BREAK 

IT  IS=0,  TEST  FOR  BREAK 

JMP 

0MPP 

NOT=0,  CONTINUE 

MP36 : 

JSR 

0 BREAK 

IS  IT  STILL  TRACKING? 

JMP 

0MPP 

.BLK 

20 

. ZREL 

; PROPAGATE  SUBROUTINE  ************* 

PROP: 

PRO 

; IT  PROPAGATES  THE  MEAN  AND 

. NREL 

; THE  VARIANCE 

PRO: 

STA 

3 ,TEM 

STORE  RTN  ADDR 

LDA 

1 , @0, 3 

LOAD  XP  OR  YP 

SUBO 

3,3 

CLR  AC0 

MOVZL# 

1,1, SNC 

TEST  SIGN 

INC 

3,3, SKP 

POS,  SET  FLAG 

NEG 

1,1 

NEG, NEGATE 

LDA 

2, A 

LOAD  A 

SUBO 

0,0 

CLR  AC0 

MUL 

PI: 

LDA 

2 , SF11 

LOAD  SF11 

DIV 

MOVZR 

3, 3, SNC 

IS  FLAG  SET? 

NEG 

1.1 

NO,  NEGATE 

LDA 

3 ,TEM 

RESTORE  RTN  ADDR 

LDA 

2,01,3 

LOAD  UX  OR  UY 

ADD 

2,1 

P2 : 

STA 

1,02,3 

STORE  XM  OR  YM 

LDA 

0,03,3 

LOAD  PP  MSW 

LDA 

1,®4,3 

LSW 

MOV 

0,0, SZR 

IS  MSW=0? 

PROl : 

JMP 

PR02-1 

NO,  JMP  PR02 

LDA 

2,B 

YES,  LOAD  B 

MUL 

LDA 

2 ,SF11 

LOAD  SF11 

DIV 

SUBO 

0,0 

CLR  AC0 

STA 

3 ,TEM 

STORE  RTN  ADDRESS 

JMP 

PR06  ; JMP  TO  ADD  G 

SUBO 

3,3 

/CLR  SHFT  COUNTER 

PR02: 

INC 

3,3 

■INC  SHFT  COUNTER 

MOVZR 

0,0, SZR  ;IS  MSW=0? 

JMP 

PRO 3 ;NO,  JMP  PR03 

MOVR 

1,1  ; YES , MOVE  LSW 

JMP 

PR04-2 

/JMP  DESHIFT 

PR03 : 

MOVR 

1.1 

■MOVE  LSW 

JMP 

PR02  ;MSW  NOT=0,  GO  BACK 

LDA 

2,B  /LOAD  B 

MUL 

o 
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PR04 : 

NEG 

3,3 

NEG  SHFT  COUNTER 

LDA 

2 , SHFT 

LOAD  SHIFT 

ADD 

2,3 

MOVZL# 

3,3, SNC 

DSHFT  LFT  OR  RGHT? 

MOV 

3,3  SNR 

NO  SHFTS? 

JMP 

PR06 

NO  SHFTS,  JMP  PR06 

JMP 

.+6 

POS,  DSHFT  RGHT 

MOVZL 

1,1 

NEG,  LFT;  MOV  LSW 

MOVL 

0,0 

MSW 

P3s 

INC 

3,3, SZR 

INC  COUNTER 

JMP 

.-3 

GO  BACK 

JMP 

PR06 

JMP  TO  PR06 

NEG 

3,3 

NEG  COUNTER 

MOVZR 

0,0 

MOVR  MSW 

MOVR 

1,1 

LSW 

INC 

3, 3, SZR 

INC  COUNTER 

JMP 

. -3 

GO  BACK 

PR06 : 

LDA 

3 ,TEM 

RESTORE  RTN  ADDR 

LDA 

2 ,G 

LOAD  G MSW 

ADD 

2,0 

ADD  MSW 

LDA 

2 ,G1 

LOAD  G LSW 

ADD 

2,1 

ADD  LSW 

STA 

0,05,3 

STORE  PXM  MSW 

STA 

1,06,3 

LSW 

JMP 

27,3 

RETURN 

.BLK 

20 

. ZREL 

; ARGUMENT  SUBROUTINE  ************** 

ARGM: 

AR 

; CALCULATES  (X**2+Y**2) /2SIGS 

.NREL 

• 

9 

AR: 

STA 

3 , TEM 

STORE  RTN  ADDR 

LDA 

1 , XM 

LOAD  XM 

MOVZL# 

1,1, SCZ 

TEST  DIGN 

NEG 

1,1 

NEG,  NEGATE 

MOV 

1,2 

DUP  XM 

SUBO 

0,0 

CLR  AC0 

MUL 

AR1 : 

STA 

0 ,TEM1 

TEMPOR  STORE  MSW 

STA 

1 ,TEM2 

AND  LSW  OF  XM**2 

LDA 

1 ,YM 

LOAD  YM 

MOVZL# 

1,1, SZC 

TEST  SIGN 

NEG 

1,1 

NEG,  NEGATE 

MOV 

1,2 

DUPLICATE  YM 

SUBO 

0,0 

CLR  AC0 

MUL 

AR2 : 

LDA 

2 ,TEM1  ; LOAD  MSN  OF  XM**2 

ADD 

2.0 

ADD 

LDA 

2 ,TEM2 

LOAD  LSW  OF  XM**2 

ADD 

2,1  ; ADD 

LDA 

2,SIGS2  ;LOAD  1*SIGS 

DIV 

J 

STA 

1 , ARG  ; STORE  EXP  ARGUMENT 

123 
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AR3 : 

JMP 

0TEM 

RETURN 

r 

. BLK 

20 

. ZREL 

;H,  HH  SUBROUTINE  **************** 

HHHH: 

HHH 

; CALCULATES  H,  AND  HH  FOR  GIVEN 

. NR£L 

; ARGUMENTS  IN  CALLING  STATEMENT 

HHH: 

STA 

3 ,TEM 

STORE  RTN  ADDR 

LDA 

1 , EXPT 

LOAD  EXPT 

LDA 

2 , IMAX 

LOAD  IMAX 

SUBO 

MUL 

0,0 

CLR  AC0 

LDA 

2 , SF11 

LOAD  SFll 

DIV 

LDA 

2 , BIAS 

LOAD  I BIAS 

ADD 

1.2 

ADD 

HI: 

STA 

2,@2,3 

STORE  HX 

LDA 

2,@0,3 

LOAD  XM  OR  YM 

SUBO 

3,3 

CLR  AC 3 

MOVZL# 

2,2, SNC 

TEST  SIGN 

INC 

3,3, SKP 

POS,  SET  FLAG 

NEG 

2,2 

NEG,  NEGATE 

SUBO 

0,0 

CLR,  AC0 

H2: 

MUL 

LDA 

2 , SIGS 

LOAD  SIGS 

DIV 

MOVZR 

3,3, SZC 

TEST  FLAG 

NEG 

1,1 

NEG,  NEGATE 

c 

LDA 

3 ,TEM 

RESTORE  RTN  ADDR 

1 

STA 

1,01,3 

STORE  HHX  OR  HHY 

i ! 

H3 : 

JMP 

13,3 

RETURN 

1 

.BLK 

20 

.ZREL 

; UPDATE  SUBROUTINE  *************; 

UPDT: 

UPD 

.NREL 

; UPDATES 

• 

t 

K AND  P 

UPD: 

STA 

3 ,TEM 

STORE  RTN  ADDR 

LDA 

0,00,3 

LOAD  PXM  OR  PYM  MSW 

LDA 

1,01,3 

LSW 

LDA 

2 , MSK2 

LOAD  SHFT  MASK 

JSR 

0SCL 

SHIFT  P 

STA 

0 , TEM13  ; STORE  SHIFTED  PM  MSW 

STA 

1,TEM14 

LSW 

STA 

3 ,TEM12 

STORE  SHFT  * 

Ul: 

LDA 

2 ,SF11 

LOAD  SFll 

DIV 

LDA 

3 , TEM  ; RESTORE  RTN  ADDR 

LDA 

2,05,3  ; LOAD  HHX  OR  HHY 

SUBO 

0,0  ;CLR  AC0 

MOVZL* 

2, 2, SNC  ;TEST  SIGN 

INC 

0,0, SKP  ;POS,  SET  FLAG 

NEG 

2,2  ;NEG,  NEGATE 

U2: 

STA 

0, SGN  ; STORE  SIGN  FLAG 

0 
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r ■ 

STA 

2 , AHH 

STORE  ABS(HHX) 

0 

SUBO 

MUL 

0,0 

■CLR  AC0 

’ 

STA 

0 ,TEM10 

STORE  H*P  MSW 

STA 

1 ,TEMll 

LSW 

LDA 

2 , SF11 

;LOAD  SF11 

U3s 

DIV 

1 

LDA 

2,  AHH 

■LOAD  ABS(HHX) 

SUBO 

MUL 

0,0 

■CLR  AC0 

1 

LDA 

2 , SF11 

■LOAD  SF11 

DIV 

LDA 

2 , TEM12 

■LOAD  SHFT  # 

U4 : 

SUBO 

0,0 

'CLR  ACO 

JSR 

@DESC 

DESHIFT  SNGL  WORD 

LDA 

2 , R 

LOAD  R 

ADD 

1,2 

•ADD;  ^DENOMINATOR 

LDA 

0,TEM10 

•LOAD  NUM  MSW 

LDA 

1,TEM11 

LSW 

DIV 

’ 

U5: 

LDA 

2 ,TEM12 

LOAD  SHFT  # 

SUBO 

0,0 

•CLR  ACO 

JSR 

0DESC 

•DESHIFT  SINGLE  WORD 

LDA 

3 ,TEM 

RESTORE  RTN  ADDR 

STA 

1 , AK 

STORE  ABS(KX) 

MOV 

1,2 

DUPLICATE  ABS(KX) 

n 

LDA 

0,SGN 

LOAD  SIGN  FLAG 

U6 : 

MOVZR# 

0,0, SNC 

TEST  SIGN 

NEG 

1,1 

NEG,  NEGATE 

STA 

1,02,3 

STORE  KX  OR  KY 

LDA 

1,  AHH 

LOAD  ABS(HHX) 

SUBO 

MUL 

0,0 

CLR  ACO 

LDA 

2 , SF11 

LOAD  SF11 

U7s 

DIV 

SUB 

1,2 

1-K*HH 

STA 

2,TEM1 

TEMP  STORE  1-K*HH 

LDA 

0 ,TEM13 

LOAD  PM  MSW 

LDA 

1,TEM14 

LSW 

LDA 

2 , SF11 

LOAD  SF11 

DIV 

08; 

LDA 

2 ,TEMl 

LOAD  1-K*HH 

SUBO 

MUL 

0,0 

CLR  ACO 

LDA 

2 ,TEM12 

LOAD  SHFT  # 

JSR 

0DESC 

DESHIFT 

LDA 

3,TEM 

RESTORE  RTN  ADDR 

STA 

0,03,3  ; STORE  PXP  MSW 

STA 

1,04,3 

LSW 

JMP 

16,3  ; RETURN 

1 0 

* 

125 

f 



. BLK 

20 

. ZREL 

;SHIFT  SUBROUTINE  ***************** 

SCL: 

SC 

; IT  SHIFTS  A DOUBLE  WORD  UNTIL  THE 

.NREL 

.MSB  IS  MASK2  (IN  AC2)  JUSTIFIED. 

SC: 

STA 

3, TAB 

STORE  RTN  ADDR 

SUBO 

3,3 

CLR  AC 3 

MOVZ 

0,0, SNR 

IS  MSW=0? 

MOVZ 

1,1, SZR 

YES,  IS  LSW*=0? 

JMP 

.+2 

ONE  IS  NOZERO 

JMP 

§TAB 

BOTH*=0 , RETURN 

MOV 

0,0, SNR 

IS  MSW=0? 

JMP 

SC  3 

YES,  JMP  SC 3 

AND# 

0,2, SNR 

SHFT  LFT  OR  RT? 

JMP 

SC2 

LEFT 

SCI: 

MOVZR 

0,0 

MOVE  RIGHT  MSW 

MOVR 

1,1 

LSW 

INC 

3,3 

INC  AC 3 

AND# 

0,2, SZR 

DONE? 

JMP 

SCI 

NO,  GO  BACK 

NEG 

3,3 

YES 

JMP 

§TAB 

RETURN 

SC2: 

MOVZR 

2,2 

SHIFT  MASK 

MOVZL 

1,1 

SHIFT  LSW 

MOVL 

0,0 

MSW 

INC 

3,3 

INC  SHFT  COUNTER 

AND# 

0,2, SNR 

DONE? 

JMP 

SC2+1 

NO,  DO  IT  AGAIN 

JMP 

0TAB 

YES,  RETURN 

SC3: 

AND# 

2,1, SZR 

16  SHFTS  OR  LESS? 

JMP 

SC2 

LESS,  JMP  SC2 

MOV 

1,0 

SHFT  LFT  BY  16 

SUBO 

1,1 

•1 

LDA 

3 ,SF7 

LOAD  SF=16 . 

JMP 

§TAB 

RETURN 

.BLK 

10 

.ZREL 

; DESHIFT 

OESC : 

DES 

; SHIFTS 

A DOUBLE  WORD  BY  # OF  BITS 

.NREL 

; SPECIFIED  IN  AC2  BY  CALL. 

DES: 

STA 

3, TAB 

STORE  RTN  ADDR 

MOV 

2, 2, SNR 

0 SHIFTS? 

JMP 

@TAB 

YES,  0 SHFTS  - RETURN 

MOVZL# 

2 , 2 ,SZC  ;LEFT  OR  RIGHT  SHIFT? 

JMP 

DES 2 ;LEFT  SHIFT 

NEG 

2,2 

RIGHT  SHIFT 

DES1 : 

MOVZR 

0,0  ; SHFT  RT  MSW 

MOVR 

1,1  ; AND  LSW 

INC 

2, 2, SZR  ; 

INC  SHFT  COUNTER 

JMP 

DES1  ;NOT  DONE,  GO  BACK 

JMP 

0TAB  ; DONE , RETURN 

DES2: 

MOVZL 

1,1 

SHFT  LFT  LSW 

MOVL 

0,0  ; 

MSW 

o 
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INC 

2,2, SZR 

DONE? 

JMP 

DES2 

‘NO,  GO  BACK  FOR  MORE 

JMP 

@TAB 

YES,  RETURN 

. BLK 

10 

. ZREL 

; CLOCK  SUBROUTINE  ************ 

CLCK: 

CLK 

; IT  SAMPLES  DATA,  UPDATES  XP, 

• NREL 

; AND  CALCULATES  RESIDUES. 

CLK: 

STA 

3 ,TEM 

STORE  RTN  ADDR 

SUBO 

0,0 

CLR  AC0 

DOAS 

0 , ADCV 

LOAD  CHNL  SELECT 

SKPBZ 

ADCV 

IS  DATA  READY? 

JMP 

.-1 

NO,  WAIT 

DOAS 

0, ADCV 

DO  IT  TO  KILL  TIME 

SKPBZ 

ADCV 

W 

JMP 

.-1 

It 

DIC 

1 ,ADCV 

LOAD  DATA 

MOV 

1,2 

DUPLICATE  DATA 

Cll: 

ADDZL 

1,1 

*4 

ADD 

1,2 

=5*DATA  (SCALING) 

STA 

2 , ZETA 

STORE  ZETA 

LDA 

1, §6,3 

LOAD  H 

SUB 

1,2 

ZETA-H=RES 

STA 

2, §0,3 

STORE  XRES 

SUBO 

0,0 

CLR  AC0 

C2 : 

MOVZL* 

2,2, SNC 

TEST  SIGN 

INC 

0,0, SKP 

POS,  SET  FLAG 

NEG 

2,2 

NEG,  NEGATE 

STA 

2, §1,3 

STORE  ABS  (XRES) 

LDA 

1,§2,3 

LOAD  KX 

MOV 

0,3 

MOV  SGN  TO  AC3 

SUBO 

0,0 

CLR  AC0 

C3s 

MOVZL* 

1,1, SNC 

TEST  SIGN 

INC 

3, 3, SKP 

POS,  SET  FLAG 

NEG 

1,1 

NEG,  NEGATE 

MUL 

LDA 

2,SF11 

LOAD  SF11 

DIV 

MOVZR 

3 , 3 , SZC 

TEST  FLAG 

C4i 

NEG 

1,1 

NEG,  NEGATE 

LDA 

3,TEM 

RESTORE  RTN  ADDR 

LDA 

2, §3, 3 

LOAD  XM 

ADD 

2,1 

STA 

1,§4,3 

STORE  XP 

LDA 

1, §5, 3 ; LOAD  SXRES 

LDA 

2 ,TRES  ; LOAD  TIME  CONSTANT 

SUBO 

0,0  ;CLR  AC0 

C5: 

MUL 

; 

LDA 

2 ,SF11  ; LOAD  SF11 

DIV 

J 

STA 

1 , TEM1  ; TEMP  STORE 

LDA 

1,31,3  ; LOAD  AX RES 
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C6: 


TEN: 

TE: 


Tl: 


OUTPT : 
OUT: 


IEST: 

IES: 


LDA 

2, GRES 

;LOAD  GAIN 

SUBO 

0,0 

CLR  AC0 

MUL 

LDA 

2 , SF11 

LOAD  SF11 

DIV 

LDA 

2 , TEM1 

RETURN  STORE 

ADD 

2,1 

STA 

1,05,3 

STORE  SXRES 

JMP 

17,3 

RETURN 

.BLK 

20 

. ZREL 

OUT 

.NREL 

STA 

LDA 

SUBO 

MOVZL# 

INC 

NEG 

LDA 

SUBO 

DIV 

MOVR 

NEG 

LDA 

LDA 

DOB 

DOA 

JMP 

• BLK 

.ZREL 

IES 


.ZREL  ; SIGN  SUBROUTINE  **************** 
TE  ; IN  1ST  10  ITERATIONS,  THE  SIGN  OF 

.NREL  ; ERROR  EST.  IS  FORCED 


STA 

3 ,TEM 

STORE  RTN  ADDR 

LDA 

1,00,3 

LOAD  AXRES 

LDA 

0, BOUND 

LOAD  BOUND 

SUBZ# 

1,0, SZC 

IS  SXRES  GT  BOUND? 

JMP 

12,3 

NO,  RETURN 

LDA 

0,81,3 

LOAD  XP 

NEG 

0,0 

NEGATE  XP 

STA 

0,01,3 

STORE  NEW  XP 

JMP 

12,3 

RETURN 

.BLK 

20 

; OUTPUT  SUBROUTINE  ************* 
; OUTPUTS  SPECIFIED  DATA  ON 
; SPECIFIED  CHANNEL 


3,TEM 

1.01.3 

3.3 

1,1, SNC 
3,3, SKP 
1,1 

2 , SFOUT 

0,0 

3, 3, SNC 
1,1 
3,TEM 
0,00,3 

0, DACV 

1, DACV 

12.3 
10 

; IMAX  EST.  SUBROUTINE  ********** 
; CALCULATES  THE  FILTERED  AVG  OF 


; STORE  RTN  ADDR 
; LOAD  DATA 
;CLR  FLAG 
;TEST  SIGN 
POS,  SET  FLAG 
NEG,  NEGATE 
LOAD  D/A  SCL  FCTR 
CLR  AC0 

IS  NEG  POS  FLAG  SET? 
NO,  NEGATE 
RESTORE  RTN  ADDR 
LOAD  CHANNEL  SELECT 
; SPECIFY  CHNL  SELECT 
; OUTPUT  D/A 
; RETURN 


.NREL 

;I  (EST)1 

-ZETA*IMAX/H 

STA 

3,TEM 

; STORE  RTN  ADDR 

LDA 

1 , ZETA 

; LOAD  ZETA 

LDA 

0,BIAS 

; LOAD  BIAS 

SUB 

0,1 

? SUBTRACT  BIAS 

LDA 

2,  IMAX 

;LOAD  IMAX 

O 
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SUBO 

0,0 

CLR  AC0 

MUL 

LDA 

2,00,3 

LOAD  H 

LDA 

3,  BIAS 

LOAD  BIAS 

SUB 

3,2 

SUBTRACT  BIAS 

DIV 

LDA 

2 ,GI 

LOAD  FLTR  GAIN 

11: 

SUBO 

0,0 

CLR  AC0 

MUL 

LDA 

2 , SF11 

LOAD  SF11 

DIV 

STA 

1 ,TEM1 

TEMP  STORE 

LDA 

1, IMAX 

LOAD  IMAX 

LDA 

2,  GIX 

LOAD  TIME  CONSTANT 

12: 

SUBO 

0,0 

CLR  AC0 

MUL 

LDA 

2 ,SF11 

LOAD  SF11 

DIV 

LDA 

2,TEM1 

RETURN  TEMP  STORE 

ADD 

2,1 

LDA 

2 ,GII 

LOAD  GII 

LDA 

0,  R 

LOAD  R 

MOVZR 

0,0 

DIVIDE  BY  128 

MOVZR 

0,0 

M 

MOVZR 

0,0 

N 

MOVZR 

0,0 

n 

MOVZR 

0,0 

n 

MOVZR 

0,0 

N 

MOVZR 

0,0 

n 

JMP 

.+1 

SAVE  ROOM  FOR  MORE 

JMP 

.+1 

DIVIDES 

ADD 

0,2 

SUBO 

0,0 

CLR  AC0 

MUL 

LDA 

2 ,SF11 

LOAD  SF11 

DIV 

STA 

1 , IMAX 

STORE  NEW  IMAX  ESTIMATE 

LDA 

3,TEM 

RESTORE  RTN  ADDR 

13: 

JMP 

1,3  ; RETURN 

. BLK 

20 

. ZREL 

REST: 

RES 

;R(EST) -FILTERED  AVG  OF  (RES**2)/9 

• NREL 

; 

RPS: 

LDA 

1,00,3  ; LOAD  AXRES 

MOV 

1,2  ; DUPLICATE 

SUBO 

0,0  ;CLR  AC0 

MUL 

1 

LDA 

2,01,3  ;LOAD  PROPRT.  CONSTANT 

DIV 

1 

LDA 

2 ,GR  ;OAD  FILTER  GAIN 

RE1: 

SUBO 

0,0  ;CLR  AC0 
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SUBZL 

1.2.SZC 

IS  I-S  NEG?  (S>I?) 

L B2: 

JMP 

eiw 

; YES  RESTART 

JMP 

0.3 

•return 

. BLK 

20 

.ZREL 

; SZETA  SUBROUTINE*********** 

SZETA: 

SZE 

.•COMPUTES  A FILTERED  AVG  OF 

.NREL 

; ZETAS  (MEASUREMENTS) 

SZE: 

LDA 

l.GSX 

LOAD  GSX 

LDA 

2.SZET 

LOAD  SZET 

SUBO 

0.0 

CLR  AC0 

MUL 

LDA 

2.SF11 

LOAD  SF11 

DIV 

STA 

1.TEM1 

TEMPOR  STORE  RSLT 

SZ1: 

LDA 

l.ZETA 

LOAD  ZETA 

LDA 

2.GS 

LOAD  GAIN 

SUBO 

0,0 

CLR  AC0 

MUL 

5 

1 

LDA 

2.SF11 

LOAD  SF11 

DIV 

LDA 

2.TEM1 

RETURN  TEMPOR  STORE 

SZ2: 

ADD 

2,1 

STA 

l.SZET 

STORE  SZET 

JMP 

0,3 

RETURN 

.BLK 

20 

c 


o 


BOUN: 


. ZREL 
BOU 
. NR£L 


; BOUND  SUBROUTINE************** 
.•COMPUTES  THE  BOUND  FOR  SIGN 
.•FORCING  DURING  CONVERGENCE 


BOU: 

LDA 

1,SLPE 

.•LOAD  SLOPE 

LDA 

2, SZET 

.•LOAD  SZET 

SUBO 

0,0 

. CLR  AC0 

MUL 

; 

LDA 

2.SF11 

; LOAD  SF11 

DIV 

• 

LDA 

2,  MIN 

.•LOAD  CONSTANT 

BOU1: 

ADD 

2,1 

• 

§ 

STA 

1, BOUND 

.•STORE  BOUND 

JMP 

0,3 

; RETURN 

.BLK 

20 

V 

.ZREL 

.•SEARCH 

SRCH: 

SRC 

; IT  USES 

1 A SPIRAL  SCAN  TO  LOCATE 

.NREL 

;THE  TARGET  BEFORE  LOCK-ON 

SRC: 

STA 

3, TAB 

; STORE  RTN  ADDR 

SR: 

LDA 

l,NO 

; INITIALIZE: 

STA 

1#  N 

; STEP  COUNTER 

SX: 

LDA 

0,N 

; LOAD  COUNTER 

STA 

0,NC 

; STORE  IN  X-COUNTER 

SX1: 

LDA 

l.STP 

; LOAN  STEP 

LDA 

2,DAX 

.•LOAD  X-OUTPUT  REGISTER 

SX2: 

ADD 

2,1 

; NEW  X-OUTPUT  REGISTER 

STA 

l.DAX 

; STORE  IT 
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JSR 

BNDRY 

DAX 

JSR 

0OUTPT 

ZER 

DAX 

.BLK 
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SX3: 

JSR 

TRSH 

DSZ 

NC 

JMP 

SX1 

ISZ 

N 

SY: 

LDA 

0 ,N 

STA 

0,NC 

SY1: 

LDA 

1,STP 

LDA 

2,  DAY 

SY2: 

ADD 

2,1 

STA 

1,DAY 

JSR 

BNDRY 

DAY 

JSR 

eOUTPT 

ONE 

DAY 

.BLK 

10 

SY3 : 

JSR 

TRSH 

DSZ 

NC 

JMP 

SY1 

ISZ 

N 

LDA 

1,STP 

NEG 

1,1 

STA 

1,STP 

JMP 

SX 

.BLK 

20 

BNDRY : 

LDA 

1, @0, 3 

MOVZL* 

1,1, SZC 

NEG 

1,1 

LDA 

2 , SFOUT 

SUBO 

0,0 

DIV 

LDA 

2 ,MSKB 

AND 

2,1, SNR 

JMP 

1,3 

SUBO 

0,0 

STA 

0,DAX 

STA 

0 , DAY 

JMP 

SR 

.BLK 

20 

TRSH: 

LDA 

0,ZETA 

MOVZR 

0,0 

SUBO 

1,1 

DOAS 

1, ADCV 

SKPBZ 

ADCV 

JMP 

.-1 

TEST  FOR  BOUNDARY 
OUTPUT  X-CHNL 


TEST  FOR  THRESHOLD 
DONE  WITH  X-AXIS? 

NO/  CONTINUE 

YES,  INC  N AND  DO  Y-AXIS 

RESTORE  NC 

ft 

LOAD  STEP  INCREMENT 
LOAD  Y-OUT  REGISTER 

STORE  NEW  Y-OUT  RGSTR 
TEST  FOR  BOUNDARY 

OUTPUT  Y-AXIS  CONTROL 


TEST  FOR  THRESHOLD 
DONE  WITH  Y-AXIS? 
NO,  CONTINUE 
INC  N 

REVERSE  STEP  DRCTN 


DO  IT  AGAIN 

LOAD  DAX  (OR  Y) 
TEST  SIGN 
NEG,  NEGATE  IT 
LOAD  D/A  SF 
CLR  ACO 

LOAD  BOUNDARY  MASK 
TEST  FOR  BOUNDARY 
NONE,  RETURN 
RESTART  X & Y AXIS 
AT  ORIGIN 


; TESTS  FOR  THRESHOLD 
; FILTER  THE  MEASUREMENTS 
;SET  A/D  CODE 
; START  CONVERSION 
; IS  IT  READY? 

;NO,  WAIT 
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DIC 

1 , ADCV 

LOAD  DATA 

MOV 

1,2 

SCALE  TO  SF11 

ADDZL 

1,1 

N 

ADD 

1,2 

ff 

MOVZR 

2,2 

DIV  BY  2 FOR  FILTERING 

ADD 

0,2 

STA 

2, ZETA 

STORE  FILTERED  ZETA 

LDA 

1 , TEST 

LOAD  THRESHOLD  BOUND 

SUBZL 

2, 1,SZC 

ZETA-TEST 

JMP 

0TAB 

ABOVE  THRSHLD,  JMP  MAIN 

JMP 

0,3 

BELOW,  RETURN 

. BLK 

20 

. ZREL 

; TABLE  STORAGE  SUBROUTINE  ****** 

TABLS: 

TABL 

; IT  STORES  1ST  264  VALUES 

.NREL 

; SPECIFIED  IN  CALLING  STMNT 

TABL: 

STA 

3 ,TEM 

STORE  RTN  ADDR 

LDA 

0 , TABC 

LOAD  TBL  COUNTER 

INC 

0,0, SZR 

IS  TBL  FULL? 

INC 

0,0, SNR 

«f 

JMP 

4,3 

YES,  RETURN 

STA 

0,TABC 

STORE  TBL  COUNTER 

LDA 

2 ,TBlX 

LOAD  X-TBL 

ADD 

0,2 

TRUE  ADDRESS 

TAB1 : 

LDA 

1,00,3 

LOAD  1ST  X VALUE 

STA 

1,0,2 

STORE  IT 

LDA 

1,01,3 

LOAD  2ND  X VALUE 

STA 

1,1,2 

STORE  IT 

LDA 

2 ,TB1Y 

LOAD  Y-TBL 

ADD 

0,2 

TRUE  TBL  ADDRESS 

LDA 

1,02,3 

LOAD  1ST  Y VALUE 

STA 

1,0,2 

STORE  IT 

LDA 

1,03,3 

LOAD  2ND  Y VALUE 

STA 

1,1,2 

STORE  IT 

TABC2 : 

JMP 

4,3 

RETURN 

TB1X: 

TBX+277 

TBlYs 

TBY+277 

.BLK 

20 

.ZREL 

; EXPONENTIAL  SUBROUTINE  ******** 

EXP: 

EX 

; USES  TABLE  LOOKUP  TO  CALCULATE 

.NREL 

; EXP  (**X) 

EX: 

STA 

3 ,TEM 

STORE  RTN  ADDR 

LDA 

2,00,3 

LOAD  ARGUMENT 

LDA 

0,MSK4 

LOAD  ARGUMENT  MASK 

AND# 

2,0, SZR 

IS  ARG>8? 

JMP 

EX3 

YES,  JMP  EX3 

ADDZL 

2,2  ; SHIFT  LEFT  BY  2 

LDA 

3 ,MSK  ; LOAD  MASK 

MOVS 

3,1  j SWAP  BYTE  INTO  AC1 

EX1: 

ANDS 

2,3  ;HOB  IN  AC3 

AND 

2,1  ;LOB  IN  AC1 

LDA 

2 ,TBL  ;LOAD  TABLE  ADDR 
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EX2 : 


EX3 : 


TBLs 

TABLE: 

3701 

3604 

3511 

3417 

3330 

3242 

3156 

3073 

3012 

2732 

2654 

2600 

2524 

2452 

2402 

2332 

2264 

2217 

2153 

2110 

2046 

2006 

1746 

1707 

1652 

1615 

1561 

1526 

1473 

1442 

1411 

1361 

1332 

1304 


ADD 

2,3 

; ACTUAL  ADDRESS 

LDA 

2,1,3 

;LOAD  ARGMNT+1 

LDA 

3,0,3 

; LOAD  ARGUMENT 

MOV 

3,0 

; DUPLICATE  VALUE 

SUB 

2,0 

; SLOPE  IN  AC0 

MOVS 

0,2 

; PREPARE  SLOPE  FOR  MUL 

SUBO 

0,0 

;CLR  AC0 

MUL 

• 

$ 

SUB 

0,3 

; EXPT  IN  AC 3 

JMP 

.+2 

• 

$ 

SUBO 

3,3 

; ARG>8 , EXPT=0 

LDA 

2 ,TEM 

; RELOAD  RTN  ADDR 

STA 

3,ei,2 

; STORE  EXPT 

JMP 

12,2 

; RETURN 

• BLK 

20 

TABLE 

4000 
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t 

1256 

213 

34 

6 

1 

( 

1231 

207 

33 

6 

1 

1204 

203 

33 

5 

1 

1161 

177 

32 

5 

1 

1135 

173 

31 

5 

1 

1113 

167 

30 

5 

1 

1071 

164 

27 

5 

1 

1047 

160 

27 

5 

1 

1026 

155 

26 

4 

1 

1006 

151 

25 

4 

1 

766 

146 

25 

4 

1 

746 

143 

24 

4 

1 

727 

140 

23 

4 

1 

711 

135 

23 

4 

1 

673 

132 

22 

4 

1 

655 

127 

22 

4 

1 

640 

125 

21 

3 

1 

623 

122 

21 

3 

1 

607 

117 

20 

3 

1 

573 

115 

20 

3 

1 

557 

113 

17 

3 

1 

544 

110 

17 

3 

1 

531 

106 

16 

3 

1 

516 

104 

16 

3 

TBX: 

. BLK 

340 

504 

102 

15 

3 

TBY: 

. BLK 

340 

472 

100 

15 

3 

.END 

0 

460 

76 

15 

3 

447 

74 

14 

2 

* 

436 

72 

14 

2 

425 

70 

13 

2 

415 

67 

13 

2 

404 

65 

13 

2 

374 

63 

12 

2 

365 

62 

12 

2 

355 

60 

12 

2 

346 

57 

11 

2 

337 

55 

11 

2 

330 

54 

11 

2 

321 

53 

11 

2 

313 

51 

10 

2 

305 

50 

10 

2 

276 

47 

10 

2 

271 

46 

10 

2 

263 

44 

7 

2 

255 

43 

7 

1 

250 

42 

7 

1 

243 

41 

7 

1 

236 

40 

7 

1 

231 

37 

6 

1 

224 

36 

6 

1 

220 

35 

6 

1 

0 

t 

I 

; 
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